Quantitative Assessment of Dynamic Stability Characteristics for Jet Transport in Sudden Plunging Motion

The main objective of this article is to present a training program of loss control prevention for the airlines to enhance aviation safety and operational efficiency. The assessments of dynamic stability characteristics based on the approaches of oscillatory motion and eigenvalue motion modes for jet transport aircraft response to sudden plunging motions are demonstrated in this article. A twin-jet transport aircraft encountering severe clear-air turbulence in transonic flight during the descending phase will be examined as the study case. The flight results in sudden plunging motions with abrupt changes in attitude and gravitational acceleration (i.e. the normal load factor). Development of the required thrust and aerodynamic models with the flight data mining and the fuzzy-logic modeling techniques will be presented. The oscillatory derivatives extracted from these aerodynamic models are then used in the study of variations in stability characteristics during the sudden plunging motion. The fuzzy-logic aerodynamic models are utilized to estimate the nonlinear unsteady aerodynamics while performing numerical integration of flight dynamic equations. The eigenvalues of all motion modes are obtained during time integration. The present quantitative assessment method is an innovation to examine possible mitigation concepts of accident prevention and promote the understanding of aerodynamic responses of the jet transport aircraft.


Introduction
The databases of flight simulators and autopilots are designed and established based on data in nominal flights in current state of art. The flight condition of transport aircraft does not include in normal design process and this situation cannot be tested easily during the aircraft certification process, but most likely to cause the aircraft loss of control in service is so-called off-nominal condition. The flight mechanism of off-nominal condition cannot be simulated by wind tunnel tests easily, and it is not possible to directly perform flight tests with real aircraft. In other words, the design value of aerodynamic performance is not consistent with the actual aerodynamic performance during the flight operations. It is appropriate that off-nominal flight conditions cannot be combined with aircraft design and flight simulators (Chang and Tan 2013). More importantly, the pilots can't have much chance to hone in this flight condition.
The adverse weather related to off-nominal flight conditions include (but are not limited to): Heavy rain, convective storms (Sauer et al. 2019), turbulence (Hamilton and Proctor 2003;Lee and Lan 2000), crosswind (Lan and Chang 2018;Weng et al. 2006), wind shear (Min et al. 2011;Wang et al. 2016), in-flight icing (Lan and Guan 2005;Thompson et al. 2004). Adverse weather impacts the flight safety of transport aircraft and operational efficiency of air traffic management. An innovation tool of commercial flight route planning program is developed to avoid the deep convective storms (Sauer et al. 2019). This program provides an overview to make a better air route based on original trajectory operations for safe, efficient, and comfortable operating flights to the airlines.
Among the different adverse weathers, clear-air turbulence is all the more important in-flight safety since it is hard to detect and predict. Clear-air turbulence is the leading cause of serious personal injuries in non-fatal accidents of commercial aircraft. One main type of motion that causes flight injuries in clear-air turbulence is the sudden plunging motion with abrupt changes in attitude and gravitational acceleration (i.e. the normal load factor). The assessment of flying qualities during sudden plunging motion for an aged twin-jet transport aircraft encountering severe clear-air turbulence through digital 6-DOF flight simulations in transonic flight was presented in reference (Chang et al. 2009). To numerically integrate the 6-DOF dynamic equations of motion and determine the eigen-modes of motion at the same time at every instant were the important process of that approach. The longitudinal equations being coupled with the lateral-directional equations and some nonlinear effects being incorporated were the advantages for 6-DOF flight simulations. The eigenvalue equations were solved with QR transformation (Chang et al. 2009). Unfortunately, that approach was difficult to identify the individual modes of motion from these eigenvalues, from one instant, to another because of the rapid changes of aerodynamic forces and moments in turbulence.
To have a better approach, this article presents flight data mining and the fuzzy-logic modeling of artificial intelligence technique to establish nonlinear and unsteady aerodynamic models for six aerodynamic coefficients based on the flight data of a new twin-jet transport aircraft. This new twinjet encountered severe clear-air turbulence twice during the descending phase. The oscillatory derivatives extracted from the aerodynamic models are then used in the study of variations in stability characteristics during the sudden plunging motion. The longitudinal and lateral-directional motion modes are analyzed through digital flight simulation based on decoupled dynamic equations of motion. The eigenvalue equations are formulated in the form of polynomials and solved. The eigenvalues for short-period, phugoid (long-period), dutch roll, spiral, and roll modes of motion are estimated based on damping ratio and undamped natural frequency. A positive real part of the eigenvalues is to indicate unstable motion of the related modes. The evaluations of unstable conditions for each eigen mode with the oscillatory derivatives are also demonstrated in the present paper.

Flight Data Mining Technology
The input data for fuzzy-logic modeling is extracted from post-flight data. There are many flight parameters recorded by QAR, but some of them are not related to the current study, such as light signs, landing gear retracting, ..., etc.. The steps are to collect and organize a large amount of data and obtain the required information from them for specific analysis and applications; this process is so called flight data mining (Han et al. 2011). The flight data mining process is divided into two parts: the developments of a non-linear and unsteady aerodynamic database and the aerodynamic models.

Development of nonlinear and unsteady aerodynamic database
This section describes the development process of the demand data selection, data frequency upscaling and data make up, compatibility check, and input information of aircraft main geometry and moment of inertia data. The flowchart of each step is shown in Fig. 2.1.

Fig. 2.1 Flowchart for development of a non-linear and unsteady aerodynamic data database
Before developing a non-linear and unsteady aerodynamic database, one must ensure that the selected parameters and corresponding data are the required data for study. The parameters and corresponding data required for current study include the following items: 1) To select the required flight status, flight control surface, and engine raw data for a specific flight number from the actual flight QAR or FDR data.
2) The main geometric and moment of inertial data include appearance size, wingspan, average chord length, and wing area, etc.. Most of the above are public data of aircraft manufacturers and can be obtained via the Internet. However, the thrust performance of the entire aircraft is non-public data and must be obtained from the Flight Planning and Performance Manual (Boeing Series) or Flight Crew Operation Manual (Airbus Series).
The frequency of each parameter sensor is different from 1Hz to 8Hz due to the properties and characteristics of the parameters. For example, attitude angles such as roll, pitch, and heading angles are 4Hz; aileron and elevator angles are 2Hz; the vertical acceleration sampling rate is 8Hz. To preserve the properties and characteristics of each parameter, it is necessary to avoid high frequency data loss, so the frequency of all original data sampling is uniformly arranged to 8 Hz.
The data frequency increase and data make up of this research is developed based on the method 5 of monotone cubic spline interpolation (Wolberg and Alfy 2002). This method is commonly used in aircraft flight test analysis to provide a smooth curve. Data upscaling and data completion of blanks are important steps for aerodynamic database to have high accuracy and reliability in applications to operational efficiency enhancement of transport aircraft.

Compatibility Check
Typically, the longitudinal, lateral, and vertical accelerations (ax, ay, az) along the (x, y, z)-body axes of aircraft, angle of attack α, Euler angles (φ, θ, and ψ), aileron deflection (δa), elevator (δe), rudder (δr), stabilizer (δs), …, etc. are available and recorded in the QAR of all transport aircraft. Since the recorded flight data may contain errors (or called biases), compatibility checks are performed to remove them by satisfying the following kinematic equations (Lan and Chang 2012;: where g is the gravitational acceleration and V is the flight speed. Let the errors be denoted by , respectively for ax, ay, az, etc. These errors are estimated by minimizing the squared sum of the differences between the two sides of Eqs. 2.1~2.6.
These equations in vector form can be written as: where super fix "-" stands for the mean value and the subscript "m" indicates the measured or recorded values. The cost function is defined as: where Q is a weighting diagonal matrix with elements being 1.0 and z   is calculated with a central difference scheme with m z  , which is the measured value of z  . The steepest descent optimization method is adopted to minimize the cost function. As the result of analysis, variables unavailable in the QAR, such as β, p, q and r, are capable of being estimated.
The above-mentioned unavailable variables in the QAR, need initial values as a basis for correction, where the angular rates, such as p, q, r, are obtained from derivatives of φ , θ , and Ψ with time by using the method of monotonic cubic interpolation. This interpolation method is used to connect the flight data into a continuous curve to obtain the slope of the curve. The value of β cannot be obtained from the derivative. The initial value of β is assumed to be 0 at the time of estimation, and it is calculated when the Eq. 2.6 is satisfied.
Aircraft main geometry and moment of inertia data The main objective of this article is to present a dynamic stability monitoring method based on the approach of oscillatory motion and eigenvalue motion modes for the twin-jet transport aircraft response to sudden plunging motions. The data of main geometry and moments of inertia are presented in Table 2.1. Moment of inertia-xz axes, Ixz ( kg•m 2 ) 0.0

Non-linear and Unsteady Aerodynamic Database
The actual physical phenomenon of an aircraft flying in the atmosphere is the physical movement along the flight trajectory over time. Since all flight variables recorded are based on the body axes, it is more convenient to estimate the force and moment coefficients for aircraft on the same axes system as follows (Roskam 2018): 1) force coefficients acting on the body axes of the aircraft (2.14) 2) moment coefficients about the body axes of the aircraft where m is the aircraft mass; q is the dynamic pressure; S is the wing reference area; Cx, Cz, and Cm are the longitudinal aerodynamic force and moment coefficients; Cy, Cl, and Cn are the lateraldirectional aerodynamic force and moment coefficients; and Ixx, Iyy, and Izz are the moments of inertia about x-, y-, and z-axes, respectively. The products of inertia, Ixy, Ixz, and Iyz, are assumed zero in the present case; but are included in the equations because non-zero values may be available in other applications. The terms, Tx, Ty, Tz, and Tm, represent the thrust contributions to the force in the direction of x-, y-, z-axes, and to the pitching moment, respectively.
In equations (2.15) ~ (2.17), pq, qr, pr are the product terms of two dependent variables or 2 2 , r p are the square terms of a variable; those are non-linear terms in mathematics. Therefore, equations (2.15) to (2.17) belong to the nonlinear differential equation system. Non-linear differential equations are very difficult to solve mathematically, and the physical phenomena of motion are more complicated (Fazelzadeh and Sadat-Hoseini 2012). In the physical sense, the product term of the two dependent variables is related to each other, which is the coupling effect of nonlinear motion. The raw flight data of QAR or FDR can be used to establish a non-linear and unsteady aerodynamic database through the steps of flowchart in Fig. 2.1.
The terms on the left-hand side of Eqs. (2.15) -(2.17) are rolling, pitching, and yawing moments, respectively. The moments due to inertia coupling on the right-hand side in Eqs. (2.15) -(2.17) will be produced due to the abrupt change in the flight attitude before and during the severe clear-air turbulence encounters.

Development of nonlinear and unsteady aerodynamic models
This section describes the development process of nonlinear and unsteady aerodynamic models, as shown in Fig. 2.2. The flowchart of Fig. 2.2 is based on the data of the non-linear and unsteady aerodynamic database; the thrust data in aerodynamic database has been separated from aerodynamic coefficienties; Model strucg-n, the numerical model structure, is setup through parameter selection, time segment selection, sampling rate setting, and then go through the refining process of fuzzy-logic modeling to establish the required models; the operational efficiency enhancement can be studied after of model defuzzification. T are the axial thrusts along body x-, y-, z-coordinate, and pitching moment about yaxis, respectively, in the equations (2.12) ~ (2.14), and (2.16). It can be known from the various thrusts that the force and moment acting on the aircraft will be affected by the engines. Both thrust and aerodynamic forces are generated together, and the two effects cannot be accurately separated only from the flight record data during analysis. To accurately estimate the aerodynamic coefficient, it is necessary to obtain an accurate thrust value to separate the thrust effects, and then the accurate aerodynamic coefficients can be obtained ).
Since the values of thrust for aircraft in flight cannot be directly measured in the current state of the art, they are neither recorded in the QAR nor FDR. The manufacturers of engines agreed that using such parameters as the Mach number, airspeed, flight altitude, temperature, rpm of the pressure compressors, and engine pressure ratios are adequate to estimate the engine thrust (Chang 2014). To study aircraft performance with fuel consumption, data from the flight manual for the fuel flow rate For the Pratt & Whitney turbofan engines, thrust (T) is defined by EPR, so that the thrust model is set up as: For GE or CFM turbofan engines, the rpm of the low-pressure compressor (N1) is used to set the level of thrust, so that the thrust model is set up as: In the present study, the twin-jet with GE turbofan engines under study will be illustrated. The actual thrust in operation is obtained by using the recorded variables in the QAR, in particular the fuel flow rates. Once the thrust models are generated as a function with the flight conditions of climbing, cruise, and descending phases; one can estimate the thrust magnitude by inserting those flight data into the aerodynamic database.
When studying the nonlinear and unsteady aerodynamic characteristics of motion, it can be found that the nonlinear and unsteady aerodynamics with hysteresis (Chang et al. 2010;Lan and Chang 2018;Paziresh et al. 2017) are affected by many parameters of motion state. The aerodynamic models developed by the fuzzy-logic modeling method is very suitable for aerodynamic research of coupling motion, and can be used to enchance the operational efficiency of transport aircraft. The development process of nonlinear and unsteady aerodynamic models are shown in Fig. 2.2.
Modeling means to establish the numerical relationship among certain variables of interest. In the fuzzy-logic model, more complete necessary influencing flight variables can be included to capture all possible effects on aircraft response to structure deformations. The longitudinal main aerodynamics is assumed to depend on the following ten flight variables (Chang et al. 2010): The The most representative of numerical differentiation theory is the difference method. This paper uses the central difference method to take the derivative. The central difference method to obtain the derivative is to take the differential point as the center point, and respectively select a value close to 0 from the left and right of the differential point to compare, and obtain the derivative of the certain point.
Assuming a function f(x), if one wants to find the derivative at the point x, take x as the center point and select the value of the distance between the left and right sides close to 0, namely f(x+) and f(x-), and then find the Taylor expansion and the principle of the central difference method are as follows: Taylor expansion to find f(x+ Taylor expansion to find f(xx ∆ )： Subtracting (2.22) from (2.23), one gets Although the above-mentioned numerical differentiation method is used to obtain the derivative of a certain point of the function curve, in fact, it can also be applied to the derivative analysis of experimental data. However, the experimental data are scattered and discontinuous points, and interpolation must be used to connect the points into a continuous curve. Fuzzy-logic modeling plays the function of interpolation.
The tangent slope of a point on the curve is the derivative value of this point. Add and subtract a small same value to the abscissa value of this point respectively, corresponding to the slope of the line connecting two points on the curve, which is calculated by the finite difference method. Fuzzylogic modeling is extremely important in the aerodynamic performance analysis of actual flight data.
Here are two examples of applying the fuzzy-logic modeling model to take derivatives: derivative of the unsteady aerodynamic model can be calculated using the central interpolation method. The formula of the central difference method is as follows: (2.25) ∆α =0.1 degrees mean that α changes up and down by 0.1 degrees, while all other variables remain unchanged.
The damping derivative (Clp) is proposed from the model Cl. The formula of the central interpolation method is as follows: The calculations of all other aerodynamic derivatives are derived from the same method in this paper.

Flight Simulation
Since the rapid changes of aerodynamic forces and moments in turbulence, it is difficult to identify the individual modes of motion from these eigenvalues through digital 6-DOF flight simulations from one instant to another. Regarding to the dynamic stability monitoring, one approach to solve this problem is to use the approximate modes of motion obtained from decoupled longitudinal and lateral-directional equations as guidance.
The decoupled linearized longitudinal equations of motion, which are decoupled from the lateral-directional motions in reference of (Roskam 2018) are given by where the u X , α X , and e X δ are the dimensional variations of force along X axis with the speed, angle of attack, and elevator angle, respectively; the other dimensional derivatives of Z and M are described and given in reference of (Roskam 2018 ; the dimensional derivatives of Y, L, and N are described and given in reference of (Roskam 2018). The characteristic equations for Eqs, (2.27) ~ (2.29), and Eqs.
(2.30) ~ (2.32) are polynomials of the 4th degree and their roots are solved with a quadratic factoring method based on the Lin-Bairstow algorithm in reference of (Hovanessian and Pipes 1969).
The 4 th degree polynomial of longitudinal characteristic equations has 4 roots; they are two complex conjugates (Lan and Chang 2012). The short-period mode is one of the two complex conjugates. Another one is the phugoid mode (long period mode). Each mode has the same real part, but with imaginary parts of equal magnitude and opposite signs. The 4 th degree polynomial of lateraldirectional characteristic equations also has 4 roots; they are one pair of complex conjugates and two real values. The pair of complex conjugates represent the dutch roll mode. One of two real values represent the spiral mode and another one represents the roll mode.

Fuzzy-logic Modeling Algorithm
Modelling procedures start from setting up numerical relations between the input (i.e. flight variables) and output (i.e. flight operations or aircraft response). To obtain continuous variations of predicted results, the present method is based on the internal functions, instead of fuzzy sets (Karasan et al. 2018), to generate the output of the model.
The fuzzy-logic algorithm can take advantage of correlating multiple parameters without assuming explicit functional relations among them (Dang et al. 2017;Shi et al. 2017). The algorithm employs many internal functions to represent the contributions of fuzzy cells (to be defined later) to the over-all prediction. The internal functions are assumed to be linear functions of input parameters as follows (Weng et al. 2006): the computed gradient tends to be small and the convergence is slow. To accelerate the convergence, the iterative formulas are modified by using the local squared errors to give: for r=0,  is weighted factor of the i th cell; and the index j of the data set, where j=1, 2, …, m, and m is the total number of the data record; and the "op" stands product operator of its elements in this paper.
Once the fuzzy-logic aerodynamic models were set up, one can input to the model influencing variables to describe the flight conditions under analysis. In the present paper, all aerodynamic derivatives are computed with central differences through the aerodynamics models.
The fuzzy-logic algorithm can take advantage of correlating multiple parameters without assuming explicit functional relations among them (Chang 2014;Chang et al. 2010). The algorithm employs many internal functions to represent the contributions of fuzzy cells (to be defined later) to the over-all prediction. The internal functions are assumed to be linear functions of input parameters as follows (Roskam 2018):

Numerical Results and Discussions
In the present study, the accuracy of the established unsteady aerodynamic models with six aerodynamic coefficients by using fuzzy-logic modeling technique is estimated by the sum of squared errors (SSE) and the square of multiple correlation coefficients (R 2 ). All the aerodynamic derivatives in the study of dynamic stability characteristics are calculated with these aerodynamic models of aerodynamic coefficients.
Cn 3 2 4 2 2 2 3 3 2 2 2 13824 0.9435 The twin-jet transport aircraft encountered clear-air turbulences in revenue flights. As a result, several passengers and cabin crews sustained injuries, this event was classified as the aviation accidents. To examine the dynamic stability characteristics, it is imperative to understand the flight environment in detail.
The corresponding flight data for twin-jet transport with plunging motion is presented in Fig. 3.2.
The dataset of time span from t=7480~7739 sec used for the modeling are extracted from the FDR.
This twin-jet transport encounters severe clear-air turbulence twice during this time span. In Fig. 3.2(a), the variations of normal acceleration (az) show the highest az being 2.05g around t = 7483 sec and the lowest being-1.05 g around t = 7484 sec in the first turbulence encounter, the highest az being 1.91g around t = 7682 sec and the lowest being -0.16 g around t = 7684 sec in the second one.  of pitch rate (q) and roll rate (p) are shown in Fig. 3.3(a); the variation of p is larger than that of q.
The yaw rate is not shown because it is small throughout. It implies that pitch angle (θ) does not vary as much as AOA and the variations of roll angle (φ) is large, especially during the second turbulence encounter.
is about a certain value with the magnitude of small fluctuation; the variations of roll angle (φ) is large; the yaw rate is not shown in Fig. 3.3(a), because it is small throughout. It implies that twin-jet transport aircraft has crosswind encounter, but this crosswind does not have big variations in directions.

Fig. 3.3 Time history of twin-jet transport aircraft response to dynamic flight motion
In general, the configuration in static stability assessment is based on enough values of steady damping derivatives to keep stable condition. The plunging motion is the biggest factor that easy causes aerodynamic instability, the configuration will produce a time-varying pressure distribution on the surface, which involves the compressibility effect. The longitudinal, roll, and directional (Cnr)osc are from positive at t= 7483.5 sec to negative at t= 7485 sec. It implies that the effects of β derivative is to cause the directional stability more unstable. In essence, the effects of α  -derivative on (Czq)osc, and -derivative on (Clp)osc are small in Figs 3.4 and 3.5. However, the effect of α  -derivative on (Cmq)osc is to improve the stability in pitch; while the effects of -derivative is to cause more directional instability. These results indicate that the turbulent crosswind has some adverse effects on directional stability and damping. Although the dynamic derivatives tend to be small for the present configuration, these are much helpful to understand the unknown factors of instability characteristics.
In the present study, the longitudinal and lateral-directional motion modes are analyzed based on the damping ratio (ζ) and undamped natural frequency (ωn). The roots of the complex conjugate are as follows: where n ζω − is real part (i.e. in-phase) and 2 1 ζ ω − ± n i are imaginary (i.e. out-of-phase) parts. λr andλi represent eigenvalues of real and imaginary parts, respectively. If λr is positive, the system is unstable; if it is negative, the system is stable (Chang et al. 2009).   For the present purpose, if the severity is defined by the lowest load factor developed, then this transport aircraft developed -1.05g with ∆h=57.3m (around t = 7483 sec). The phugoid mode in Fig.   3.6(b) has all positive real eigenvalues; it implies insufficient in longitudinal oscillatory damping and with virtual mass effects during sudden plunging motion. So as this transport does not have sufficient lateral-directional oscillatory damping. Therefore, it may be concluded that the plunging motion is more severe because of its unstable condition.

Concluding Remarks
The main objective of this paper was to present a dynamic stability monitoring method based on the approach of oscillatory motion and eigenvalue motion modes for jet transport aircraft response to sudden plunging motions. Three obvious achievements could be concluded from the present paper as follows: 1) The flight data mining and the fuzzy-logic modeling techniques were shown to be effective in establishing the nonlinear and dynamic aerodynamic models through the flight data of FDR. The fuzzy-logic aerodynamic models had robustness and nonlinear interpolation capability to generate the required independent variables in predicting the oscillatory motion during severe clear-air (c) turbulence encounter.
2) The eigenvalues of both longitudinal and lateral-directional motion modes were analyzed through digital flight simulation based on decoupled dynamic equations of motion. The decoupled concept was an innovative consideration for the flight simulation of six-DOF.
3) The results of flight simulation showed that the Phugoid, Dutch roll, and spiral modes were in unstable conditions for the present study transport aircraft. Those unstable conditions were easily judged by the positive real part of the eigenvalues during sudden plunging motion.
The oscillatory motion and eigenvalue approaches for the dynamic stability assessments were much helpful to understand the effects of severe clear-air turbulence encounter on dynamic stability characteristics. The simulated results could provide the training program of loss control prevention for the airlines to enhance aviation safety and operational efficiency.

Acknowledgments
This work was supported in part by the Joint Funds project of the National Natural Science φ, θ, ψ = Euler angles in roll, pitch, and yaw, respectively, deg.