Abstract

Examining the flight quality of an aircraft to ensure the stability of the aircraft, increase maneuverability, and make the aircraft easier to control by the pilot necessitates an examination of the natural stability of the system. Within the scope of the paper, the frequency domain response of the F-16 aircraft dynamics is analyzed using Simulink models considering two different flight regimes because the frequency-domain methods have many distinct and important advantages over time-domain methods. Aerodynamic, propulsive, and atmospheric databases are used to create the nonlinear model. The trim analysis for cruise flights is carried out to obtain trim parameters. The aircraft is numerically linearized using the small perturbation theory. The linearized dynamics for each trim condition are used to create transfer functions for each input. The linear model is subsequently examined in the frequency domain to obtain information about the dynamic behavior of the aircraft, and flight quality analysis was examined by considering the lateral and longitudinal modes of the aircraft by international standards. It has been clearly understood the stability augmentation system design has critical importance for the modes with unstable or long steady-state duration.

Keywords

• Frequency Domain Analysis
• F-16 Aircraft
• Lateral Motion
• Longitudinal Motion
• System Dynamics

1. Introduction

Stability analysis is a critical phenomenon that needs to be addressed for achieving the targeted mission according to the aircraft type. Passenger comfort, the pilot’s ability to control the aircraft, the calculation of flight performance, how accurately the sensors will work, and many other criteria (Hess, 2007) can be revealed by analyzing the natural stability of the aircraft. Especially with the developing technology, enhancing technical capacities and travel time maneuverability, and increasing time-varying data such as ammunition and fuel necessitates performance and efficiency calculations and high-accuracy tests (Dündar et al., 2020). These tests are carried out using the aerodynamic efficiency data arising from the design of the aircraft, the equations of motion calculated from the aerodynamic coefficients and the effectiveness of the control surfaces, flight performance calculations, etc.

Frequency domain analysis is useful for dynamic analysis, robust controller design (Hess, 2007) and verification, guidance and trajectory studies, air combat research, and many other missions (Garza and Morelli, 2003). It indicates how well the aircraft responds to a range of input frequencies. Changes in amplitude and phase fit into a specific pattern that depends on the aircraft’s design, from which information about the system is revealed. Also, the higher bandwidth means the faster-commanded speed or position adjustment time (Wescott, 2006). The bode is an analysis plot for a linear time-invariant system in the frequency domain that includes a magnitude and a phase plot (Atangana and Akgül, 2020).

Morelli et al. mentioned that (Morelli and Grauer, 2020; Morelli and Jared, 2015) frequency-domain studies for F-16 fighter aircraft have crucial advantages over time-domain studies due to the reasons such as noise rejection, accurate modeling, the modeling of parametric uncertainty, enabling simple and effective approach, data extraction, and validation, real-time modeling, flying qualities modeling, and computational efficiency. Along with these advantages, a disadvantage of frequency-domain estimation is associated with the practical limitation of identification techniques to linear dynamical systems only (Klein, 1978 and 1980). Fu stated (Fu, 1990) that the frequency response of the transfer function facilitates the enhancement of uncertain transfer matrices for multi-input, multi-output systems (MIMO). Meanwhile, these results are utilized to determine the H∞ norm, gain margin, and phase margin and to improve the diagonal dominance of the MIMO indeterminate system. A frequency-domain study is also used to determine the stability parameters of the aircraft (Morelli, 2000, Millidere, 2021, Shukla et al., 2017). The short-period characteristic evaluated in the longitudinal motion was obtained by driving the elevator input at different frequencies. With some estimated aerodynamic derivatives of the system, the second order transfer function and accordingly the damping ratio and natural frequency data can be obtained (Milliken, 1947, Greenberg, 1951).

In this paper, the frequency domain response of the F-16 fighter airplane dynamics is analyzed with performed Simulink models considering altitude of 3000 feet (914.4 m) and 30000 feet (9144 m) trimming conditions. To acquire comprehensive nonlinear dynamics, the plant model of the F-16 fighter airplane is established by modeling aerodynamic, gyroscopic, gravitational, and propulsive force and moments concerning atmospheric effects. The modeling part comprises various data such as the geometric configuration of fighter airplanes, aircraft wind tunnel test aerodynamic and propulsion data, given inertia-mass parameters, and engine angular momentum.

In kinematic calculations, Euler angular representation is chosen to obtain transfer functions and state-space representations to make the final results understandable. The most important advantage of quaternion is that it eliminates the singularity problem. However, it is more difficult to interpret quaternion than Euler angles. Since there is no singularity at the points where linearization is made, it does not affect the analysis made in the frequency domain, and in order to make the final results understandable, both transfer functions and state-space representations are made according to Euler angles.

The aircraft operation intervals, which also define where the trim conditions are, should be specified prior to linearization. The intervals can be defined for take-off, cruise, or another mode. The states and required control inputs can then be derived by trim analysis of these operating points. Once the trim variables are obtained, the nonlinear dynamics of the aircraft can be linearized to obtain the state-space equations. In this paper, reappraisal is carried out in this study by looking at the dynamics in the frequency domain. In the end, the aircraft dynamics are also investigated purpose whether they are appropriate for human and aircraft dynamic response interfaces, as successful maneuvering, approaching, and tracking missions are all dependent on them. As a result, for the given criteria, this analysis should also be carried out in terms of the damping ratio, time constant, and natural frequency of lateral and longitudinal modes.

This paper is organized as follows. The nonlinear equation of the system dynamics is addressed in the second chapter. Then, the trim conditions for cruise flight and linearization process of the nonlinear model are expressed to obtain transfer functions and break frequencies which constitutive they are tabulated. Also, the state-space model of the F-16 aircraft is composed. In the third chapter, Bode diagrams and pole maps are demonstrated. Finally, obtained results and their effects on stability modes are analyzed and interpreted.

2. Equation of Motion

Fig. 1.    Earth-fixed and body-fixed reference frames (Stevens and Lewis, 1992)

Considering the equation of motion of an aircraft, various assumptions are made in order to simplify the equations. First, the airframe is assumed to be a rigid body, so the center of gravity of the aircraft won’t have a relative velocity for the Body-fixed reference frame. Second, the mass of aircraft is assumed to be constant, and finally, mass distribution remains constant during the flight, resulting in a constant inertia tensor. The equation of motion of an aircraft is mainly obtained on Earth-fixed and Body-fixed reference frames. These reference frames are represented in Fig. 1.

Six nonlinear equations of motion were utilized to model nonlinear aircraft dynamics for translational and rotational motion (Stevens and Lewis, 1992). These equations are shown in (1-4).

(1)

(2)

(3)

External forces acting on the aircraft in a body-fixed reference frame, such as aerodynamic, gravitational, and propulsive forces, are represented by , , and . Sometimes, these equations are used by transforming them into earth-fixed or wind axes reference frames. For the earth-fixed frame, transformation can be done by using the relation below.

(4)

Where ϕ, θ, and ψ are the Euler angles that correspond to orientation in the earth-fixed frame and the transformation matrix , provides a transformation from body-fixed frame to earth-fixed frame. The other transformation to wind axes can be done using the (5).

(5)

On the other hand, rotational equation of motion can be found using Newton’s second law,

 (6)

Hence, angular velocity rates can be found as

(7)

(8)

(9)

In general applications, regarding aircraft geometry, some product inertia terms that are  and  may be neglected. Another transformation can be performed for rotational kinematics, which incorporates angular body rates and Euler angle rates,

(10)

3. Modeling

During the flight, some external forces and moments are acting on the aircraft. There are variety of factors that cause these forces and moments. These variables and model types will be discussed in this section. Aerodynamic, propulsion, atmospheric, gyroscopic, and gravitational models are all covered in this article.

For aerodynamic modeling, Nguyen et al. (1979) provided extensive data about the fighter aircraft F-16, and the data that is used in this paper is obtained from that article by Nguyen (1979). This data consists of two main parts, which are breakpoints and data. The data part is based on the data which is obtained from wind-tunnel testing of F-16 (Nguyen, 1979). On the other hand, breakpoints store the angle of attack, slip angle, or surface deflection data in which wind-tunnel testing is conducted.

(11)

(12)

Where the dynamic pressure is q ̅, the wing reference area is S, the span length is b, and the aerodynamic chord length is C.

For propulsion modeling, In F-16 aircraft, propulsion is provided by a turbofan jet engine. The engine data is given in reference to Nguyen (1979) for F-16 aircraft, and it is used to calculate thrust force . Regarding provided data, thrust force depends on Mach number, altitude of the aircraft, and throttle settings of idle, mil, and max. Thrust value is calculated by interpolating provided data with input variables. An interpolation must be done considering power percentage to find thrust value within these points. If the throttle level is within the 50-0 percent,

(13)

where ts is the percent throttle setting within 0-1, and if the throttle level is within 50-100, thrust becomes,

(14)

propulsive forces can be defined as,

(15)

the moment can be calculated using relation below,

(16)

 represents the angle between the z-axis of reference frames and  represents the angle between the y-axis of propulsion and body reference frame. are assumed to be zero, which makes the propulsive force vector  has an only component in the x-axis of the body frame.

Rotating rigid bodies have angular momentum. Suppose that an external moment acts on the rigid body that generates angular velocity. In that case, gyroscopic moments are generated to conserve previous angular momentum by producing the counter moment called the gyroscopic precession effect. According to the reference Nguyen (1979), the angular momentum of the aircraft is considered as 216.9 kg.m²/sec in the x-direction of the propulsion system reference frame. Therefore, in this system, gyroscopic moments can be found as,

(17)

Due to gravitation, there is always applied gravitation force on aircraft. As the attitude of the aircraft changes, the gravitational vector in the body axis changes. This can be modeled by using the earth-to-body reference frame transformation matrix  and gravity vector in the earth-fixed reference frame. Force components can be found as

(18)

For actuator modeling, some time lag can be observed when a command was given to an actuator. As a result, these actuator dynamics are also incorporated into the simulation system to make it more realistic. The actuator models are considered first order system. Rate and position limiters are also added to the model. Additionally, the references of actuator models can be given as Nguyen, 1979 and Sonneveldt 2006.

Fig. 2.    General Simulink model of the F-16 aircraft before trimming phase

4. Trim Conditions and Linearization

Linearization of equations of motion and decomposition of linearized equations into two sets, longitudinal and lateral, are studied in this section. The operational points that determine the numerical outcome of the trim conditions must be determined before linearization. Trim conditions will be found by finding the point where the forces and moments acting on aircraft body are zero using the optimization algorithm. The cost function to be minimized is specified in (19). The constraint equations that will set the limits of optimization are translational dynamics, rotational dynamics, translational kinematics, rotational kinematic, and some additional physical definitions for specific of the trim conditions. Translational dynamics are introduced in (1-3), translational kinematics are introduced in (4), rotational dynamics are introduced in (7-9), and rotational kinematics are introduced in (10).

(19)

States are V_T, β, α, (u, v, w), p, q, r, φ, θ, ψ, h, γ and control Inputs are δ_e, δ_a, δ_r, δ_T . Four states are defined in question V_T (velocity), h (altitude), γ (gamma angle), and φ angle. The trim is mostly searched for the defined condition. ζ_specified = [ V_T, h, γ, φ]^T, V_T = constant_, h = constant, γ = constant, φ = 0,  ζ_unknown=[ β, α, p, q, r, θ, ψ, δ_e, δ_a, δ_r, δ_T]^T. Thus, the cost function J(ζ_specified, ζ_unknown) for the general straight flight can be solved with 11 unknowns­­ for 11 equations. The constraint equations are translational dynamics, rotational dynamics, translational kinematics, rotational kinematic, and some additional physical definitions for specific trim conditions. Simulation of F-16 fighter aircraft is performed regarding these trim analyses. The states and control inputs are initialized considering findings in trim analysis, and the dynamics of the aircraft are analyzed.

Optimization algorithms make minimization according to the targeted subject. According to the selected optimization algorithm, local minimum or global minimum can be found. The cost function is primarily used to evaluate the current performance of the model by comparing the predicted results. For the J cost function to be equal to zero, all the functions that make it up must be zero. While the selected optimization algorithm solves this, it tries to bring the cost function closer to zero. Scatter-Search-Based global optimization was used in MATLAB Global Search Method to find trimmed states and inputs. A Global Search object presents a solver that attempts to locate the solution with the lowest objective function value. Global Search solver first generates trial points employing a Scatter Search method. These trial points are then filtered and fmincon is started from each of the filtered points. For case 1, the fmincon function converges to 1.535628e-11; furthermore, the Global Search algorithm converges to 5.6463e-14, corresponding to the least minimum value. For case 2, fmincon function converges to 0.00602 where corresponds to the local minimum. Since the value found with fmincon is not the global minimum, a different algorithm has been tried. The minimum cost function 3.5680e-13 was found via the Global Search algorithm. The selected algorithm cannot converge any further. While searching the global minimum, the following constraints and initial values were referred.

Table 2.   Trimmed state derivatives (x)  of Case 1

Table 4.   Trimmed Inputs (u) for Case 1

Table 5.   Trimmed state derivatives (x) of Case 2

Table 7.   Trimmed Inputs (u) for Case 2

The six nonlinear differential equations that govern aircraft motion are nonlinear differential equations. That can be solved with a variety of numerical integration techniques. However, it is practically difficult to obtain closed solutions (equations for each variable). Closed solutions to the dynamic behavior of the aircraft can provide helpful information (linear models). Hence, alternative strategies such as obtaining partial derivatives of the equation of motion or small perturbation approaches (Jafarov and Tasaltm, 1999) are used to linearize the aircraft’s equations of motion. By applying the conditions in Table 2 and Table 3, the

Table 1.    Constraints and initial values of Case-1 and Case-2

Table 3.   Trimmed States (x) of Case 1

Table 6.   Trimmed States (x) of Case 2

Linear state-space model of longitudinal dynamics for case 1 in (20)

(20)

Linear state-space model of lateral dynamics for case 1 in (21).

(21)

Linear state-space model of longitudinal dynamics for case 2 in (22).

(22)

Linear state-space model of lateral dynamics for case 2 in (23).

(23)

Table 8.   Transfer functions of trimmed points

Table 9.   Break frequencies of transfer functions

linearization process is completed, and these equations of motion are mathematically separated into two sets as longitudinal and lateral motion. After linearization, the equations of motion are divided into longitudinal motion, including variables 𝑉𝑡, 𝛼, 𝑞, 𝜃, and lateral motion, including states 𝛽, 𝜙, 𝜓, 𝑝, 𝑟. With thanks to the perturbed linearization algorithm, state-space models of Case-1 and Case-2 were generated with 2¹¹ and 2¹⁰ precision, respectively. The relevant models can be expressed as follows.

The known flight conditions were chosen with respect to constant altitude and wings level as α= q, β= -y, and angular position rates p = q = r = 0.  Trimmed states and inputs are calculated in Table 2 and Table 3.

Frequency domain analysis is done using Bode diagram for single input single output systems. For this reason, longitudinal and latitudinal dynamics expressed in state-space form should be converted into transfer functions. Table 4 shows the F-16 aircraft’s transfer functions. Table 5 also shows the break frequencies of transfer functions.

5. Frequency Domain and Flying Qualities Analysis

Transfer functions and root locus plots are provided for various modes. The outcomes of the open-loop study may be classified to represent various aspects of flight handling performance. Two flight conditions will be examined based on their behavior in different modes. Stability augmentation system and control augmentation system designs are required based on the data inspected from the frequency domain response. The design is expected to meet the desired dynamics, and these dynamics are determined by both military and civil aviation standards. All these design and performance requirements are specified in MIL-F-8785C (Flying Qualities of Piloted Airplanes), which contains these military standards (MIL-F-8785C, 1980).

The criteria in the standards have varying constraints depending on the aircraft’s class, flight phase, and flying quality. Therefore, it is necessary to determine the class of the F-16 fighter aircraft. F-16 aircraft can be classified as Class IV, Category A, or Level 1 according to MIL-F-8785C specifications. The limits that the aircraft dynamics are determined by looking at this classification (Atak, 2020). The damping ratio of the F-16 should be between 0.35< <1.30, and its natural frequency should be between 0.28 < < 3.6. The damping ratio of the short-term mode should be as great as possible. If there is insufficient damping, serious problems may emerge. For the phugoid mode, the damping ratio  ≥ 0.04 and   requirements are specified in the MIL-F-8785C standards. The frequency domain and flight quality analyses are shown in the bode diagrams below for the two situations studied.

Fig.3.     Bode diagrams for Case-1 and Case-2.

Fig. 4.    Bode diagrams for Case-1 and Case-2.

The frequency response of wind velocity V_T to horizontal tail input for case-1 and case-2 is shown. There is a pole for case-1 and zero for case-2 on the right side of the imaginary axis. Such systems are known as non-minimum phase systems. While the system in case-1 is unstable, the system in case-2 is stable. For both systems, at low frequencies, the gain margin value and phase shift remain constant. In case 2, a small peak in gain is observed when the input frequency crosses the short period break frequency, and the gain continues to decrease gradually at higher frequencies. At higher frequencies, the phase lag related to the short period mode results in a constant total phase lag of -180 degrees. The response steadily decreases as the control input frequency is raised, with an increasing phase lag. In particular, the dominant effect of Case-2 belongs to the phugoid mode, where the gain peaks because of the low damping ratio. It is explicitly located inside the observable bandwidth.

The frequency responses of the angle of attack  to horizontal tail input  for case-1 and case-2 are shown. The system in case-1 is unstable because of its positive pole.  can not be seen because of the workspace of the figure. Because of the positive pole, case-1 is an unstable and nonminimum phase system, whereas case-2 is a stable and minimum phase system. The gain margin value is constant for both cases at extremely low frequencies. When the input frequency reaches the phugoid break frequency, the case-2 gain plot shows a peak that implies a small damping ratio, and the gain continues to decrease slowly at higher frequencies. There is significant phase lag induced by phugoid and phase lead caused by  between the frequencies of 0.07 and 0.2. As the control input frequency is increased, the response diminishes steadily with increasing phase lag. In both cases, the system parameters are derived using the transfer function for pitch angle, where the input is a horizontal tail.

Fig. 5.   Bode diagrams for Case-1 and Case-2.

Fig. 6.    Bode diagrams for Case-1 and Case-2.

The frequency response of side pitch angle , to horizontal tail input  for case-1 and case-2 are shown. For case-1, one pole is placed on the right side of the imaginary axis. Thus, this system is an unstable and nonminimum phase system. For case-2, there is not any zero or pole on the right side of the imaginary axis. Thus, this system is a stable and minimum phase system. The gain margin value remains constant at very low frequencies. As the control input frequency is increased, the attitude response decreases continuously with increasing phase lag after break frequency  which caused a small peak on the gain margin. The case-1 results of these 3 transfer functions show that the longitudinal motion of the F 16 aircraft is unstable. But the longitudinal motion of the F 16 aircraft is stable for case-2.

The frequency response of sideslip angle , to aileron input  for case-1 and case-2 are shown. There is no pole on the right half side of the root locus diagram. But it has three positive zeros. Thus, this system is a nonminimum phase system. The right half plane zero has gained similar to that of the left half plane zero, but its phase nature is like a pole i.e., it adds adds negative phase to the system. Instead, the phase increases from 0 to 90 degrees, its phase increases from 0 to -90 degrees. This causes a delay in a system response which can lead to instability if not addressed. There is a peak at the dutch-roll frequency. But the gain is reduced by around −5 dB, which means that the pilot would see no significant oscillatory sideslip behavior. For case 2, there is no pole on the right half side of the root locus diagram. But it has three positive zeros. Thus, this system is a nonminimum phase system. Complex zeros at 0.30 rad/s frequency contribute approximately 15 dB into gain margin.

Fig. 7.     Bode diagrams for case-1 and case-2.

Fig. 8.     Bode Diagrams for Case-1 and Case-2.

The frequency response of roll rate p, to rudder input  for case-1 and case-2, are shown. Two systems are nonminimum phase systems because of positive zero values. Despite these zero values, systems are stable. The dutch-roll resonant peak in gain and subsequent roll-off in both gain and phase in case-1 are both standard and readily explained. If the effects at extremely low frequencies are neglected, it is not unrealistic to suggest that the usable bandwidth is a bit greater than the dutch-roll mode frequency. For case-2, increasing input frequency is affected obviously at dutch-roll break frequency. Dutch-roll frequency also adds around 10 dB to the system. In comparison between case-1 and case-2, case-1 has a less dominant spiral mode pole than case-2. As a result, the gain of case-1 ramps up till the spiral mode frequency, while the other’s gain ramps down.

The frequency response of yaw rate r, to rudder input  for case-1 and case-2 are shown. These 2 systems are minimum phase systems because there is not any positive value on the right side of the imaginary axis. In case 1, the phase lag increases in spiral mode dynamics at extremely low input frequencies until the effect of the second order numerator becomes apparent. The gain plot shows a steady but significant decrease with increasing frequency to reach a minimum of −50 dB at 𝜔_n𝜓, the resonant frequency of the second order numerator factor. The gain rises rapidly with a further increase in frequency to reach a maximum of -25 dB at the dutch-roll frequency. For case-2, the gain plot shows a steady but significant decrease with increasing frequency to reach a minimum of −50 dB at 1/T_r roll mode break frequency. The gain rises rapidly with a further increase in frequency to reach a maximum of -31.2 dB at the dutch-roll frequency. C-1 and C-2 refer to Case-1 and Case-2 in Table 10.

Fig. 9.   Bode Diagrams for Case-1 and Case-2.

Table 10. All frequency domain values of longitudinal and lateral dynamics of F16 aircraft

Mode analysis can be used to examine the aircraft’s longitudinal and lateral transfer functions in addition to the frequency response of the systems. The derived transfer functions explain how the system responds to control input or disturbances, and the characteristic equation of transfer functions defines the characteristics of the system. These characteristic equations are fourth-order polynomial equations that each root corresponds to a different mode in aircraft analysis.

The transfer function poles correspond to phugoid or short period dynamic stability modes in longitudinal dynamics. The short term dynamic mode is a damped oscillation about y-axis on the body-fixed reference frame. The modes have a relatively high damping ratio and damped frequency compared to phugoid mode. Changes in velocity in the x-direction on a Body-fixed reference frame are not greatly affected since it is described as an oscillation in the y-axis. However, primary characteristics that are affected dominantly by the short-period mode are body angular rate in x direction, pitch angle and the velocity in the z-axis on body-fixed reference frame. The phugoid mode, on the other hand, may be represented as an oscillation in velocity u. The mode has a low damping ratio as well as a low damped frequency. So that, compared to the short-period mode, its poles are closer to the imaginary axis, and it has a slower response. Additionally, in this mode, velocity changes are coupled with pitch angle and altitude changes. Furthermore, with small fluctuations, the angle of attack tends to remain consistent. On the x-axis, the phugoid mode dominates the states, pitch angle, and velocity.

Roll mode, spiral mode, and dutch-roll mode are the three dynamic stability modes in lateral dynamics. Roll mode or subsidence roll can be defined as rolling dynamics around the x- axis in the body reference frame. The roll mode is a non-oscillatory mode that has a single pole in the real axis on a complex plane. It also affects the roll angle and roll rate. Another dynamic mode is spiral mode, also defined as non-oscillatory dynamics in lateral direction. When the sideslip angle disturbance excites the aircraft, the aircraft fin affects the roll angle, and dihedral effect influences the roll angle. It can be affected by disturbances on the sideslip angle. Therefore, regarding this, it can be stated that is a coupled motion (Bajodah at al., 2018) in yaw, sideslip, and roll. On longitudinal dynamics, the dutch-roll mode is identical to the short-period mode, but it occurs in the z-axis on the body-fixed reference frame. It usually has complex conjugate poles on the complex plane. This dynamics couples with the roll angle and the sideslip angle.

The poles in the longitudinal study are -0.124±0.0715j, -1.72, and 0.0571. Real poles belong to the phugoid mode. One of the poles is positive, which means it is unstable. Because there is instability, in this case, the case does not meet any criteria. For the second case, the poles of phugoid mode are at -0.00201±0.0827j with 0.0243 damping ratio and 0.0827 rad/s natural frequency. The damping ratio is less than 0.04; therefore, Level 1 criteria are not met. However, the criteria for Level 2 is a damping ratio greater than 0; therefore, it is satisfied. This means although the response is relatively slow, the oscillatory response in phugoid mode eventually stabilizes.

The poles for the short period mode in the first flight condition are -0.124±0.0715j for short period frequency and dampening. The short period mode has a damping ratio of 0.867 and a natural frequency of 0.143 rad/s. The natural frequency of 0.143 does not meet any Level 1, 2, or 3 criteria for category A flights. However, it satisfies the requirements for some factor sensitivity parameters in category B. The damping ratio satisfies the criteria of being at least 0.05 for flight phase category B. The short-period modes for the second flight condition are in the -0.413±0.364j, with a damping ratio of 0.75 and a natural frequency of 0.551 rad/s. The damping ratio is sufficient once more, and the natural frequency is greater than in the previous case. The behavior satisfies the requirements of flight phase category B for all aircraft classes and C for aircraft class II-L and III for level 2 regarding standard MIL-STD-1797A. However, it does not fit categorie A and C for aircraft I, II-C, and IV for level 2.

For spiral mode, minimum time to double amplitude condition of lateral mode, the pole for the spiral mode is located in -0.0185, which is very close to the imaginary axis that makes this mode very slow to develop in the real case. Any real value that is smaller than 0.693 / T_2min is in compliance with spiral mode requirements. Where T_2min is the minimum time to double amplitude. Since the pole here is negative, it is guaranteed to be smaller than the required value. For the second flight condition, the scenario is similar. Although the pole magnitude is smaller, -0.006, it still lies on the left-hand side of the required value. This means that all the levels of spiral mode handling are achieved.

For maximum roll mode time constant, the pole of the roll mode is -3.19 for the first flight condition. The compliance region for roll mode root is . The  takes the values of 1.0, 1.4, 3.0, and 10 for different handling quality levels, yet -3.19 remains on the compliance region for any of them, so the behavior meets all the level’s criteria. However, the second flight condition’s roll mode pole is -0.848, which means while still meeting most of the requirement T_rmax values, it does not satisfy the requirements when the T_rmax is 1.0. Since the root would lie in the right-hand side of this value, outside of the compliance region. Therefore, in this condition, the handling qualities meet the criteria for Level 2 and Level 3. However, there are some shortfalls for Level 1.

For the minimum dutch-roll frequency and damping, the complex conjugate poles of dutch-roll modes are located in  with 0.139 damping ratio and 2.75 rad/s natural frequency. The damping ratio falls within satisfactory region for the Level 3 and Level 2 qualities; however, it only meets the requirement for B and C flight categories for Level 1. The minimum natural frequency is quite enough for each level, but the damping ratio is the limiting factor here. For the second condition, the complex conjugate poles of the dutch-roll mode are . The damping ratio of the dutch-roll mode is 0.119, and the natural frequency is 2.05 rad/s. Again, the natural frequency is quite enough, while the damping ratio still falls short of meeting the Level 1 category A criteria. Handling qualities are classified as the same for both conditions.

6. Conclusions

A Simulink model is developed in this paper to simulate F-16 dynamics. The aircraft’s equation of motion is also simulated. Some properties like Mach number, altitude change with aircraft dynamics, and aircraft’s actuator dynamics are also included in the system designed to make it more realistic. Then, the nonlinear aircraft dynamics are linearized around the trim conditions. The trim conditions of the aircraft are tried to be found as much as accurately using global search algorithms. Because it is also important for the accuracy of the linearized model. After finding the linearized model, the transfer functions are analyzed for frequency response and flying quality criteria. Along with the analysis in the frequency domain, mode analysis has a vital role, especially for flight performance. The phugoid, short period, roll mode, spiral mode, and dutch-roll modes were investigated longitudinally and latitudinally, respectively. It has been seen that the stability augmentation system and control augmentation system design has a critical importance for the modes with unstable or long steady-state duration. In order to make these designs, classification, level, and category of the aircraft must be determined in accordance with MIL-STD-1797A standards. For case 1, it is found that the longitudinal transfer function is unstable, which is caused by aircraft parameters coupling with different atmospheric conditions. On the other hand, for the second case, both lateral and longitudinal linearized dynamics are found as stable condition.

The results obtained for the flight at these two different altitudes reveal the effect of altitude on flight performance. It is also used as performance evaluation criteria in analyzes such as fuel consumption, controller performance, and actuator saturation of the change in flight performance. The stability augmentation system and autopilot design for F-16 aircraft will be the topic of further research after this study.

CRediT Author Statement

Abdurrahim Bilal Özcan: Conceptualization, Methodology, Software, Data curation, Writing- Original draft preparation. Elbrus Caferov: Supervision, Validation, Writing- Reviewing and Editing.

Nomenclature

References

Atangana, A., Akgül A., 2020. Can Transfer Function and Bode Diagram Be Obtained from Sumudu Transform. Alexandria Engineering Journal, 59(4), pp.1971-1984  https://doi.org/10.1016/j.aej.2019.12.028.

Atak, K.U. 2020, Aircraft Model and Flight Control System Design, Master Thesis, Department of Control and Automation Engineering, Istanbul Technical University, İstanbul, June.

Bajodah, A. H., Mibar, H., Ansari, U. 2018, Aircraft Motion Decoupling of Roll and Yaw Dynamics Using Generalized Dynamic Inversion Control. 26th Mediterranean Conference on Control and Automation (MED), Zadar, Croatia, pp.1–9. doi:  https://10.1109/MED.2018.8442505.   

Dündar, Ö., Bilici M., Ünler T., 2020. Design and performance analyses of a fixed wing battery VTOL UAV,

Engineering Science and Technology, an International Journal, Volume 23, Issue 5, pp. 1182-1193, DOI: https://doi.org/10.1016/j.jestch.2020.02.002.

Fu M. 1990. Computing the frequency response of linear systems with parametric perturbation. Systems & Control Letters, 15(1), pp.45-52, https://doi.org/10.1016/0167-6911(90)90043-T.

Garza, F. R., Morelli E. A., 2003, A Collection of Nonlinear Aircraft Simulations in MATLAB. NASA Langley Research Center, Hampton, Virginia, available at: https://www.cs.odu.edu/~mln/ltrs-pdfs/NASA-2003-tm212145.pdf (accessed 22 March 2022).

Greenberg, H., 1951, A Survey of Methods for Determining Stability Parameters of an Airplane from Dynamic Flight Measurement. Ames Aeronautical Laboratory Moffett Field, Calif, Washington, available at: https://ntrs.nasa.gov/api/citations/19930082979/downloads/19930082979.pdf (accessed 22 March 2022).      

Hess, R.A. 2007, Frequency-Domain Design/Analysis of Robust Flight Control Systems. In System Control Technologies, Design Considerations & Integrated Optimization Factors for Distributed Nano UAV Applications,  Available at: https://www.sto.nato.int/publications/STO%20Educational%20Notes/RTO-EN-SCI-175/EN-SCI-175-03.pdf (accessed 22 March 2022).

Jafarov, E.M., Tasaltm, R., 1999. Design of Longitudinal Variable Structure Flight Control System for the F-18 Aircraft Model with Parameter Perturbations, Kohala Coast, USA , (), pp.607–612. https://doi:10.1109/cacsd.1999.808716

Klein V., 1978, Aircraft Parameter Estimation in Frequency Domain. American Institute of Aeronautics and Astronautics 4th Atmospheric Flight Mechanics Conference, https://10.2514/6.1978-1344.

Klein V., 1980, Maximum Likelihood Method for Estimating Airplane Stability and Control Parameters from Flight Data in Frequency Domain. NASA Technical Paper 1637. The George Washington University, Langley Research Center, available at: https://ntrs.nasa.gov/api/citations/19800015831/downloads/19800015831.pdf (accessed 22 March 2022).

MIL-F-8785C, 1980, Military Specification-Flying Qualities of Piloted Airplanes.

Millidere M., 2021, Optimal Input Design and System Identification for an Agile Aircraft. A Thesis Submitted to The Graduate School of Natural and Applied Sciences, Middle East Technical University.

Milliken, W. F., 1947. Progress in Dynamic Stability and Control Research, Journal of the Aero. Sci., Vol. 14, pp.494 519. 

Morelli, E. A., 2000. Real-Time Parameter Estimation in the Frequency Domain. Journal of Guidance, Control, and Dynamics, 23(5), pp.812–818. https://doi:10.2514/2.4642.

Morelli E. A., Grauer J. A., 2020. Practical Aspects of Frequency-Domain Approaches for Aircraft System Identification, Journal of Aircraft 21(3), pp.268–291.

Morelli E. A., Jared C., 2015. Frequency-Domain Method for Automated Simulation Updates Based on Flight Data. Journal of Aircraft, 1(14). https://10.2514/1.C033121.

Nguyen, L. T., Marilyn E., Ogburn M.E., Gilbert W.P., Kibler K.S., Brown P.W., Deal P.L., 1979, Simulator study of stall/post-stall characteristics of a fighter airplane with relaxed longitudinal static stability (Vol. 12854). National Aeronautics and Space Administration, available at: https://ntrs.nasa.gov/api/citations/19800005879/downloads/19800005879.pdf, (accessed 22 March 2022).

Shukla, V. V., Nusrath K., 2017. Modeling & Identification of High Performance Aircraft in Frequency Domain. Journal of Mechanical and Aeronautical Engineering Research, 1(1), pp 38 – 45.

Sonneveldt L., 2006, Nonlinear F-16 model description, Control &Simulation Division, Faculty of Aerospace Engineering, Delft University of Technology.

Stevens, B.L. and Lewis, F.L., 1992. Aircraft Control and Simulation. John Wiley and Sons Ltd. ISBN 0471613975.

Stevens, B. L., Lewis, F. L., Johnson, E. N., 2015. Aircraft Control and Simulation: Dynamics, Controls Design, and Autonomous Systems, 3rd Edition. John Wiley and Sons. Hoboken, New Jersey.

Wescott T., 2006. Applied Control Theory for Embedded Systems, A volume in Embedded Technology, DOI: https://doi.org/10.1016/B978-0-7506-7839-1.X5000-4.