Nomenclature
- $A\!\left( {{x_r}} \right)$
-
membership function for input variable ${x_r}$
- a z
-
vertical acceleration, g
- b
-
wing span, m
- C x , C z , C m
-
longitudinal aerodynamic force and moment coefficient
- C y , C l , C n
-
lateral-directional aerodynamic force and moment coefficients
- $\bar c$
-
mean aerodynamic chord, m
- g
-
gravity acceleration, m/s2
- h
-
altitude, m
- I xx , I yy , I zz
-
moments of inertia about x-, y-, and z-axes, respectively, kg·m2
- I xy , I xz , I yz
-
products of inertia, kg·m2
- k 1 , k 2
-
longitudinal and lateral-directional reduced frequencies, respectively
- L, M, N
-
moments acting about the (x,y,z)-body axes of aircraft, respectively, N·m
- M
-
Mach number
- m
-
aircraft mass, kg
- ${\dot m_f}$
-
fuel flow rate, kg/h
- p, q, r
-
body-axis roll rate, pitch rate, and yaw rate, deg/s
- PMic
-
pitching moment due to inertia coupling, N·m
- RMic
-
rolling moment due to inertia coupling, N·m
- YMic
-
yawing moment due to inertia coupling, N·m
- $\bar q$
-
dynamic pressure, kpa
- R 2
-
square of multiple correlation coefficients
- S
-
wing reference area, m2
- T x
-
thrust term along the x-body axis of the aircraft, N
- T m
-
thrust-moment term in the pitching equation of motion, N
- t
-
time, s
- T, W
-
thrust and aircraft weight in flight, N, respectively
- X, Y, Z
-
forces acting on the aircraft body-fixed axes along x-, y- and z-axes, N, respectively
- $\alpha $ , $\dot \alpha $
-
angle-of-attack, deg and time rate of change of angle-of-attack, deg/s, respectively
- $\beta $ , $\dot \beta $
-
sideslip angle, deg and time rate of change of sideslip angle, deg/s, respectively
- $\delta_a $ , $\delta_e $ , $\delta_r $
-
control deflection angles of aileron, elevator and rudder, respectively, deg
- $\phi $ , $\theta $ , $\psi $
-
Euler angles in roll, pitch and yaw, respectively, deg
1.0 Introduction
Clear-air turbulence has long been a difficult issue for the aviation community because it is not only the leading cause of serious injury in non-fatal accidents, but also most likely to cause the aircraft loss of control. Historically, commercial aircraft pilots encountering en-route turbulence have been using the voice communication channel to provide in-situ warnings in verbal reports to other nearby aircraft and to air traffic control [Reference Chang, Ye, Lan and Guan1]. Air travel is considered by many to be safer than riding in a car or a train, but that doesn’t mean airline passengers don’t sustain injuries – passengers sustain an inflight injury aboard commercial flights far more often than one might expect. The Federal Aviation Administration (FAA) estimates that roughly 58 people sustain inflight injuries due to clear-air turbulence every year [Reference Hedlund2]. Other people are injured from improperly stowed baggage falling from overhead bins. In nonfatal accidents, inflight turbulence is the leading cause of injuries to airline flight attendants and passengers. Experience indicates that most flight injuries in clear-air turbulence have been caused by sudden plunging motion with a localised region of strong turbulence.
According to IATA annual report in 2005, for more than $\$ $ 65.8 million USD/year lost to airlines of turbulence-related injuries to cabin crew [3]. In recent years, the transport aircraft response to severe clear-air turbulence has been of great concern among all kinds of hazardous weathers. The flight safety prevention in avoidance of injuries for passengers and cabin crews is essential for the airlines. The influence of clear-air turbulence results from abrupt changes of flight attitude and potentially can cause the transport aircraft out of control [Reference Tan, Chang, Lan and Guan4]. These effects are quite difficult to be simulated by wind-tunnel tests or conventional flight simulators.
IATA – Loss of Control In-flight (LOC-I) program [5] has made a great contribution to arrange international training courses for the pilots in recent years. IATA’s mission is to represent, lead and serve the airline industry. Membership of IATA amounts to some 290 airlines in 120 countries. The attendees of the LOC-I training course are the representatives of the Crew Evaluation Section, Standardisation Department, Flight Operation Division from airlines. The Standardisation Department of the airlines will extend the training course in their own country, if it is necessary.
In numerical meteorology, wind with horizontal dimensions ranging from around five to several hundred kilometers can be resolved adequately. However, to accurately diagnose aircraft vertical load acceleration due to the impact of vertical wind shear in clear-air turbulence, resolved scales of motion as small as 50m are needed [Reference Proctor, Hamilton and Bowles6]. In fact, spatial resolution should be well below the aircraft length of which the order is 30m [Reference Buck, Bowles and Newman7]. Such vertical wind speed of high frequency (low wavelength) is not measurable quantitatively with the existing onboard weather radars. Sheu and Lan [Reference Sheu and Lan8] estimated the unknown vertical wind speed of clear-air turbulence on a response aircraft based on the raw flight data of flight data recorder (FDR) with accuracy consistent with the measurement of the alpha sensor. The results were used to explain the associated physical phenomena in the plunging motion. Unfortunately, a lack of the measurement data of the vertical wind speed sensor on board verifies the estimated values of damping term. Otherwise, it may be possible to define the severity of the plunging motion based on the vertical plunging equation of motion.
Regarding the vertical wind impact, the following is a review of literature worth mentioning in the study of vertical load acceleration or vertical-force in-situ turbulence. In reference to Cornman et al. [Reference Cornman, Morse and Cunning9], the clear-air turbulence severity was estimated in real time from in-situ aircraft measurements. The purpose was to calculate the eddy dissipation rate from the measurement of aircraft vertical acceleration based on the assumption of a homogeneous, continuous, linear von Kármán turbulence model coupled with a linear transfer function of the aircraft dynamics. Stewart’s research [Reference Stewart10] ignored the measurement of eddy current dissipation rate, to use the vertical-force in-situ turbulence algorithm and directly reporting (or deriving) the aircraft response (hazard) instead, but it still based on steady linear aerodynamics. Proctor & Hamilton [Reference Proctor, Hamilton and Bowles6] believed that the damage of turbulence to aircraft mainly comes from the vertical wind impact, not the eddy current dissipation rate, radar reflection ability or other variables. The method of their turbulent motion prediction is based on von Kármán’s dynamic spectral density; this method is unrealistic in most applications, since real turbulence is not isotropic.
The attendees of the LOC-I training course are the airline pilots. They care about how to mitigate the hazards and formulate preventive actions based on the existing system of their transport aircraft, but the vertical wind field on a response aircraft is unknown during the clear-air turbulence encounter. How to formulate preventive actions is the problem for the pilots. A new study of movement mechanisms for a plunging motion based on the flight data of accident events in FDR is undertaken. The comparative analyses of flight environment and aircraft response to severe clear-air turbulence for two four-jet aircraft are employed. The one with higher dropped-off altitude during the plunging motion between the two is chosen to construct the movement mechanism. The non-linear unsteady aerodynamic model of the chosen transport aircraft is established through flight data mining and fuzzy-logic modeling techniques based on the flight data of the FDR (flight data recorder). The influences of varying vertical wind and crosswind on loss of control, and sources of significant increase in angles of attack will be presented in this paper. To formulate preventive actions, the situation awareness of varying crosswind encounter for the operational pilot before and during severe atmospheric turbulence encounter will be the further study task in the future. The present study is initiated to examine possible mitigation concepts of accident prevention for the purpose of providing a valuable lecture of international training courses for IATA – Loss of Control In-flight (LOC-I) program.
2.0 Numerical method development
The input data for fuzzy-logic modeling is extracted from post-flight data. There are many flight variables recorded by the FDR, but some variables are not related to the current study, such as light signs, landing gear retracting and so on. The steps are to collect and organise a large amount of data and obtain the required information from them for specific analysis and applications; this process is called flight data mining [Reference Han, Kamber and Pei11]. The flight data mining process is divided into three parts: the developments of a non-linear unsteady aerodynamic database, platform development of root cause summary in aviation accident (including aerodynamic-model setup) and fuzzy-logic modeling algorithm.
2.1 Development of nonlinear unsteady aerodynamic database
This section describes the development process of data manipulation, compatibility check, input information of aircraft main geometry and moment of inertia data, equations of motion and unsteady thrust model. The flowchart for development of a non-linear and unsteady aerodynamic database is shown in Fig. 1. The wind field is an additional section due to the study of transport aircraft response to the impact of turbulent air zone in the present paper.
2.1.1 Data manipulation
Before developing a non-linear and unsteady aerodynamic database, one must ensure that the selected variables and corresponding data are the required data for study. The variables and corresponding data required for the current study include the following items:
-
(1) To select the required raw data of flight status, wind field, flight control surface and engine for a specific flight number from the actual flight Quick Access Recorder (QAR) or FDR data.
-
(2) The main geometric and moment of inertial data include appearance size, wingspan, average chord length, wing area and so on. Most of the above are public data of aircraft manufacturers and can be obtained via the Internet or flight operation manual. 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).
2.1.2 Wind field
Wind speed affects weather forecasting, operations of aircraft and maritime and countless other implications. Wind speed, or wind flow velocity, is a fundamental atmospheric rate. Wind speed is caused by air moving from high pressure to low pressure, usually due to changes in temperature. Wind direction in the QAR or FDR is based on cardinal direction of navigation. It refers to the horizontal angle between the clockwise direction and the target direction from the true north direction of a certain point. The angle between heading and track is known as the drift angle [12].
Wind speed, wind direction and drift angle are the required variables of wind field in actual flight data of QAR or FDR. The required variables names, referred units, and referred sampling rates are shown in Table 1. These three variables are low frequency in 1Hz.
When winds are not parallel to or directly with/against the line of travel, the wind is said to have a crosswind component; that is, the force can be separated into two vector components — a crosswind component and a headwind or tailwind component [13]. The crosswind component is computed by multiplying the wind speed by the sine of drift angle. For example, a 10-knot wind coming at 45degrees of drift angle will have a crosswind component of 10knots × sin (45°) or about 7.07knots. The headwind component is computed in the same manner, using cosine instead of sine. Wind speed, wind direction and drift angle are available on the Primary Flight Display (PFD). If the magnitude of drift angles is large and with variations, the pilots should be educated that the varying crosswind encounter will occur.
A tailwind is a wind that blows in the direction of travel of an aircraft, while a headwind blows against the direction of travel. A tailwind increases the speed of an aircraft and reduces the time required to reach its destination, while a headwind has the opposite effect. Tailwind and headwind are commonly estimated in relation to the speed of an aircraft.
Ground speed can be determined by the vector sum of the aircraft’s true airspeed and the current wind speed with direction; a headwind subtracts from the ground speed, while a tailwind adds to it. Winds at other angles to the heading will have components of either headwind or tailwind as well as a crosswind component.
A crosswind is any wind that has a perpendicular component to the line or direction of travel. They affect the aerodynamics of many forms of transport. In aviation, a varying crosswind is the component of wind that blows across the runway, making landings and take-offs more difficult than it blows straight down the runway. If a varying crosswind is strong enough, it may cause runway veer-off event after touchdown.
The transport aircraft in today’s world does not have the sensor to measure the vertical wind speed. Sheu and Lan [Reference Sheu and Lan8] estimated this value on a response aircraft based on the raw flight data in QAR or FDR with the measurement of the alpha sensor as follows.
where the subscript m stands for motion and w for wind speed (including vertical wind speed).
From another point of view, an aircraft itself is in fact a big sensor in the atmospheric environment. Penetrating a turbulent air zone, the fast-varying variables of aircraft responds in a definite way depending on the imposed wind field and aircraft aerodynamic characteristics. Typically, the longitudinal, lateral and vertical accelerations (a x , a y , a z ), Euler angles ( $\phi $ , $\theta $ , and $\psi $ ) and angle-of-attack α are influenced by the wind field. In other word, the raw flight data of fast-varying variables recorded in QAR or FDR includes the effects of the wind.
Usually, a x , a y , a z along the (x, y, z)-body axes of aircraft, α, $\phi $ , $\theta $ , and $\psi $ , aileron deflection ( $\delta_{\text{a}} $ ), elevator ( $\delta_{\text{e}} $ ), rudder ( $\delta_{\text{r}} $ ), stabiliser ( $\delta_{\text{s}} $ ), and so on are available and recorded in the QAR or FDR 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 [Reference Lan, Wu and Yu14, Reference Chang and Lan15]:
where g is the gravitational acceleration and V is the flight speed. Let the errors be denoted by ${b_{{a_x}}}$ , ${b_{{a_y}}}$ , ${b_{{a_z}}}$ , ${b_p}$ , ${b_q}$ , ${b_r}$ , ${b_V}$ , ${b_\alpha }$ , ${b_\beta }$ , ${b_\theta }$ , ${b_\varphi }$ , ${b_\psi }$ , respectively for a x , a y , a z , etc. These errors are estimated by minimising the squared sum of the differences between the two sides of Equations (2)–(7).
These equations in vector form can be written as:
where overbar “–” 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 $\dot{\mathop z\limits^{\rightharpoonup}} $ is calculated with a central difference scheme with ${\mathop z\limits^{\rightharpoonup}}_m$ , which is the measured value of ${\mathop z\limits^{\rightharpoonup}}$ . The steepest descent optimisation method is adopted to minimise the cost function. As the result of analysis, variables unavailable in the QAR or FDR, such as $\beta $ , p, q and r, are capable of being estimated.
The above-mentioned unavailable variables in the QAR or FDR need initial values as a basis for correction, where the angular rates such as p, q, r are obtained from derivatives of $\varphi $ , $\theta $ and $\psi $ 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 $\beta $ cannot be obtained from the derivative. The initial value of $\beta $ is assumed to be 0 at the time of estimation, and it is calculated when the Equation (7) is satisfied.
2.1.3 Aircraft main geometry and moment of inertia data
The main objective of this paper is to examine the possible mitigation concepts of accident prevention, the nonlinear flight controllability models based on the post-flight data are established through flight data mining and the fuzzy-logic modeling of artificial intelligence technique. The data of main geometry and moments of inertia for this transport aircraft are presented in Table 2.
2.1.4 Equations of motion
The actual physical phenomenon of an aircraft flying in the atmosphere is the 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. The following equations of motion are derived based on the textbook of Roskam [Reference Roskam16]:
-
(1) force coefficients acting on the body axes of the aircraft
(13) \begin{equation}m{a_x} = {C_x}\bar qS + {T_x}\end{equation}(14) \begin{equation}\;m{a_y} = {C_y}\bar qS + {T_y}\end{equation}(15) \begin{equation}m{a_z} = {C_z}\bar qS + {T_z}\end{equation} -
(2) moment coefficients about the body axes of the aircraft
(16) \begin{equation}{C_l}\bar qSb = {I_{xx}}\;\dot p - {I_{xz}}\dot r + qr\!\left( {{I_{zz}} - {I_{yy}}} \right) - {I_{xz}}pq\end{equation}(17) \begin{equation}{C_m}\bar qS\bar c = {I_{yy}}\dot q + rp\!\left( {{I_{xx}} - {I_{zz}}} \right) + {I_{xz}}\!\left( {{p^2} - {r^2}} \right) - {T_m}\end{equation}(18) \begin{equation}{C_n}\bar qSb = - {I_{xz}}\;\dot p + {I_{zz}}\dot r + pq\!\left( {{I_{yy}} - {I_{xx}}} \right) + {I_{xz}}qr\end{equation}
where m is the aircraft mass; $ \bar q$ is the dynamic pressure; S is the wing reference area; b is the wing span; $\bar c$ is the mean aerodynamic chord; C x , C z and C m are the longitudinal aerodynamic force and moment coefficients; C y , C l and C n are the lateral-directional aerodynamic force and moment coefficients; and I xx , I yy and I zz are the moments of inertia about x-, y- and z-axes, respectively. The products of inertia, I xy , I xz and I yz , 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 T x , T y , T z and T m represent the thrust contributions to the force in the direction of x-, y-, z-axes, and to the pitching moment, respectively. T m is a thrust-moment term about y-axis, so the sign convention is negative in Equation (17)
In Equations (16)–(18), pq, qr, pr are the product terms of two dependent variables or ${p^2},{r^2}$ are the square terms of a variable; those are non-linear terms in mathematics. Therefore, Equations (16)–(18) 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 [Reference Lan, Wu and Yu14]. 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 are time-dependent recorded data; it can be used to establish a non-linear and unsteady aerodynamic database through the steps of Sections 2.1.1 to 2.1.4.
The terms on the left-hand side of Equations (16)–(18) are rolling, pitching, and yawing moments, respectively. The moments due to inertia coupling on the right-hand side in Equations (16)–(18) will be produced due to the abrupt change in the flight attitude before and during the severe clear-air turbulence encounters.
The third term on the right-hand side of Equation (16) is the rolling moment (RMic) due to inertia coupling as follows:
The product of qr in Equation (19) represents pitch rate coupling with yaw rate, so as, the pitching moment (PM ic) and yawing moment (YM ic) due to inertia coupling are also presented as follows:
2.1.5 Unsteady thrust model
${T_x}$ , $\;{T_y}$ , $\;{T_z}$ , ${T_m}$ are the axial thrusts along body x-, y-, z-coordinate, and pitching moment about y-axis, respectively, in the Equations (13)–(15), and (17). 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 [Reference Lan, Wu and Yu14].
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. A realistic thrust model is quite complex and cannot be represented by any simple equation. Since such thrust model is not available for the present study, a numerical model ties up the recorded engine data and flight operation manual to develop the thrust model by using the fuzzy-logic algorithm.
The manufacturers of engines agreed that using variables such as the Mach number, airspeed, flight altitude, temperature, rpm of the pressure compressors and engine pressure ratios are adequate to estimate the engine thrust [1]. To study aircraft performance with fuel consumption, data from the flight manual for the fuel flow rate ( ${\dot m_f}$ ) at various altitude (h), weight (W), Mach number (M), calibrated airspeed (CAS), engine pressure ratio (EPR), in cruise flight are utilised. The flight operation manual has the required lists and graphs of engine efficiency changes with altitude and speed; these data can be used to complete the engine’s thrust change curves.
For the Pratt & Whitney (PW) turbofan engines, thrust (T) is defined by EPR, so that the thrust model is set up as:
For General Electrical (GE) or CFM International (CFM) turbofan engines, the rpm of the low-pressure compressor (N 1 ) is used to set the level of thrust, so that the thrust model is set up as:
In the present study, Aircraft A with GE turbofan engines and Aircraft B with PW one, Aircraft B will be illustrated. The actual thrust in operation is obtained by using the recorded variables from engine data in FDR. The engine’s thrust change curves tabularised in numerical datasets to associate with the recorded data in FDR are the input data files for thrust model establishment.
Theoretically, clear-air turbulence (i.e. random change in u, w (or α) and v (or $\beta $ )) affect the engine performance through its effects on static and dynamic distortions at the engine face. However, its effects are not known and cannot be estimated, and therefore are ignored in the present application.
Once the thrust model is generated as a function of h, W, M, CAS, EPR, and ( ${\dot m_f}$ ) ith the flight conditions of climbing, cruise, and descent, one can estimate the thrust magnitude by inserting those flight data into the aerodynamic database.
2.2 Platform development of root cause summary in aviation accident
This section describes the development of root caused analysis in aviation accident, as shown in Fig. 2.
The flowchart of Fig. 2 is based on the data of the non-linear unsteady aerodynamic database. The fuzzy-logic modeling technique is used to set up aerodynamic models. The non-linear unsteady aerodynamic model can provide the stability and controllability derivatives to be used for the event analyses. The derivative indices can help to pinpoint the major cause more efficiently while proceeding event or accident investigation, to judge about how difficult it was for the pilot or the autopilot system to control the aircraft in severe clear-air turbulence.
2.2.1 Non-linear unsteady aerodynamic models
When studying the non-linear unsteady aerodynamic characteristics of motion, it can be found that the non-linear unsteady aerodynamics with hysteresis [Reference Schwithal, Rohlf, Looye and Liersch17, Reference Pfnur and Breitsamter18, Reference Min, Lee, Cho, Jo, Yang and Lee19] are affected by many variables of motion state. The aerodynamic models developed by the fuzzy-logic modeling method is suitable for aerodynamic research of coupling motion, and can be used to enchance the operational efficiency of transport aircraft.
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 [Reference Lan, Wu and Yu14]:
The coefficients on the left-hand side of Equation (24) represent the coefficients of axial force (C x ), normal force (C z ) and pitching moment (C m ), respectively. The variables on the right-hand side of Equation (24) denote the angle-of-attack (α), time rate of angle-of-attack (dα/dt, $\dot \alpha $ ), pitch rate (q), longitudinal reduced frequency (k 1 ), sideslip angle ( $\beta $ ), control deflection angle of elevator ( $\delta_e $ ), Mach number (M), roll rate (p), stabiliser angle ( $\delta_s $ ) and dynamic pressure ( $\bar q $ ). These variables are called the influencing variables. The roll rate is included here because it is known that an aircraft encountering hazardous weather tend to develop rolling which may affect longitudinal stability. The variable of dynamic pressure is for estimation of the significance in structural deformation effects.
For the lateral-directional aerodynamics [Reference Lan, Wu and Yu14],
The coefficients on the left-hand side of Equation (25) represent the coefficients of side force (C y ), rolling moment (C l ) and yawing moment (C n ), respectively. The variables on the right-hand side of Equation (25) denote the angle-of-attack (α), sideslip angle ( $\beta $ ), roll angle ( $\phi $ ), roll rate (p), yaw rate (r), lateral-directional reduced frequency (k 2 ), control deflection angle of aileron ( $\delta_a $ ), control deflection angle of rudder ( $\delta_r $ ), Mach number (M), the time rate of angle-of-attack ( $\dot \alpha $ ) and the time rate of sideslip angle ( $\dot \beta $ ).
2.2.2 Aerodynamic and damping derivatives
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 centre 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 centre 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 + Δx):
Taylor expansion to find f(x−Δx):
Then,
Subtracting (26) from (28), one gets
Therefore,
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. Fuzzy-logic modeling is extremely important in the aerodynamic performance analysis of actual flight data. The following are two examples of applying the fuzzy-logic modeling model to take derivatives:
The longitudinal stability derivative (C mα ) comes from the model of C m . The aerodynamic derivative of the unsteady aerodynamic model can be calculated using the central interpolation method. The formula of the central difference method is as follows:
Δα = 0.1 degrees mean that α changes up and down by 0.1 degrees, while all other variables remain unchanged.
The damping derivative (C lp) is proposed from the model C l. The formula of the central interpolation method is as follows:
where Δp is in unit of deg/s.
The calculations of all other aerodynamic derivatives are derived with the same method in this paper.
The effectiveness criteria for derivative values of stability, damping and controllability are shown in Table 3. The root caused analysis for aircraft accident due to the pilots in loss of control are based on assessment of stability, damping and controllability in flight dynamics. The magnitude values of these three derivatives are the judgement criteria for effectiveness.
Remarks: The effectiveness is judged based on the positive or negative of derivative value.
2.3 Fuzzy-logic modeling algorithm
Since the aerodynamic models are established by using flight data, modeling technique is important and need to be carefully considered. Factors that affect the modeling procedures include the mathematical tool to set up a system model and the method to identify variables of a model structure. Modeling procedures start from separating the input data into many groups, and nonlinear relations are set up between each input-output data space. To obtain continuous variations of predicted results, the present method is based on the internal functions, instead of fuzzy sets [Reference Dang, Vermeiren, Dequidt and Dambrine20], to generate the output of the model.
The internal functions are defined to cover the ranges of the influencing variables (i.e. input variables). The ranges of the input variables are all transformed into the domain of [0, 1]. The membership grading also ranges from 0 to 1.0, with “0” meaning no effect from the corresponding internal function, and “1" meaning a full effect. These internal functions are assumed to be linear functions of input variables as follows [Reference Shi, Wang, Yang, Jian and Bi21, Reference Lan and Chang22]:
Where $p_r^i$ , r = 0, 1, 2,…, k, are the coefficients of internal functions ${y_i}$ , and k is number of input variables. In Equation (33), ${y_i}$ is the estimated aerodynamic coefficient of force or moment, and X r are the variables of the input data.
The recorded data in QAR or FDR, such as flight altitude (h), calibrated airspeed (CAS), angle-of-attack (α), accelerometer readings (a x , a y , and a z ) and Euler angles ( $\theta $ , $\phi $ , and $\psi $ ), etc. is chosen as the variables to form the data for specific fuzzy models. In the present section, ${y_i}$ is defined to be an estimated aerodynamic coefficient of force or moment, and X r are the variables of the input data. The numbers of the internal functions (i.e. cell numbers) are quantified by the total number of membership functions.
The values of each fuzzy variable, such as the angle-of-attack, are divided into several ranges, each of which represents a membership function with $A\!\left( {{x_r}} \right)$ as its membership grade. One membership function from each variable constitutes a fuzzy cell. For the i th cell, the corresponding membership grades are represented by $A_r^i\!\left( {{x_r}} \right)$ , r = 1, 2,…, k. In other words, the membership functions allow the membership grades of the internal functions for a given set of input variables to be assigned.
The flowchart of refining process for fuzzy-logic modeling is presented in Fig. 3. This flowchart includes input, normalisation, membership function calculation, internal function calculation, etc. The model refining process of modeling can be treated as a system. For a given system with input variables $\;\left( {{x_1},{x_2}, \cdots\! ,{x_r}, \cdots {x_k}} \right)$ of one data point, the recorded values of each input variables are normalised by using
Hereafter ${x_{r,norm}}$ is denoted by x r for simplicity in description. The range, (x r,max − x r,min ), represents the scaling factor and usually is assumed to have a larger range than what actually appears in the data with numerous data points to be more generally applicable for the resulting model.
The membership functions partition the input space into many fuzzy subspaces, which are called the fuzzy cells. The total number of fuzzy cells is $n = {N_1} \times {N_2} \times \cdots \times {N_r} \times \cdots \times {N_k}$ . For a variable X r , the number of membership function is ${N_r}$ . Each fuzzy cell is in a different combination from others formed by taking one membership function from each input variable.
In each fuzzy cell, the contribution to the outcome (i.e. the cell output) is based on the internal function, Equation (33). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm. The output is estimated by the centre of gravity method. For the j th input ( ${x_{1,j}},\;{x_{2,j}}, \cdots\!,{x_{r,j}}, \cdots\!,{x_{k,j}}$ ), the output is as follows:
In Equation (34), $op\!\left[ {A_1^i\!\left( {{x_{1,j}}} \right), \cdots\!,A_k^i\!\left( {{x_{k,j}}} \right)} \right]$ is the 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 records. The symbol “op” stands for product operator of its elements in the present paper.
There are two main tasks involved in the fuzzy-logic modeling process. One is the determination of coefficients of the linear internal functions. The other is to identify the best structure of fuzzy cells of the model, i.e. to determine the best number of membership functions for each fuzzy variable. The coefficients are calculated with the gradient-descent method by minimising the sum of squared errors (SSE) [Reference Pfnur and Breitsamter18]:
On the other hand, the structure of fuzzy cells is optimised by maximising the multiple correlation coefficients (R 2 ):
where ${\hat y_j}$ is the output of the fuzzy-logic model, y j is the measured data, $\bar y$ is the average value of extracted dataset, and m is the total number of data points.
The aerodynamic or flight control model is defined by the values of p r i –coefficients. These coefficients are determined by minimising SSE (Equation (35)) with respect to these coefficients. Minimisation is achieved by the gradient-descent method with an iteration formula defined by:
where ${\alpha _r}$ is convergence factor or step size in the gradient method; subscript index t denotes the iterative sequence.
After simplification, Equation (37) becomes
Since the computed gradient tends to be small with Equation (38) and involves matrix iteration so that the convergence is slow and time-intensive. To accelerate the convergence, the iterative formulas are modified by using the local squared errors to give:
For r = 0,
and for r = 1, …, k,
In Equations (39) and (40), , $op[A_1^i\!\left( {\left( {{x_{1,j}}} \right), \cdots\!,A_k^i\!\left( {{x_{k,j}}} \right)} \right]$ 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.
Equations (39) and (40) are based on a point-iteration in variables identification algorithm. The iteration during the search sequence stops, when one of the following three criteria [Reference Pfnur and Breitsamter18] is satisfied:
-
(1) Cost = SSE t = minimised
-
(2) $RER = \dfrac{{SS{E_t} - SS{E_{t - 1}}}}{{SS{E_t}}} \lt {\varepsilon _2}$
-
(3) t = t max
In the above criteria, SSE t is the sum of squared errors (SSE) being the cost function and RER = (cost_current-cost_previous)/cost_current to be denoted by “RER” (i.e. the relative error) for simplicity in descriptions; ${\varepsilon _2}$ is the required precision criteria. The algorithm variables identification plays the role of refining model. The point-iteration in variables identification algorithm is also shown in the flowchart of Fig. 3.
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.
3.0 Numerical results and discussions
Two transport aircraft under the study encountered severe clear-air turbulence in revenue flights. As a result, several passengers and cabin crews sustained injuries, because of which these two events were classified as the aviation accidents. To investigate the root causes of the aviation accidents, it is imperative to understand the aerodynamic and flight environments first, and then study the aircraft responses in operations. The flight dynamic mechanism that induced the sudden plunging motion can be constructed from the study of crosswind influence, and sources of significant increase in angles of attack.
3.1 Aerodynamic and flight environments
The variations of main flight variables in aerodynamic and flight environments are presented in Table 4. The time spans of severe clear-air turbulence encountered are underneath each aircraft as shown in Table 4. The corresponding flight data of the Aircraft A with severe clear-air turbulence encountered at the altitude around 11,277m in transonic flight. The variation of vertical acceleration, the highest a z is 1.67g and the lowest is 0g. The variation of α is approximately in phase with a z . The highest α is about 6.2deg and the lowest is −5.9deg. At the same time, M is dropped from around 0.88 to 0.83. The variations of wind speed is from +79.8 to 40.9m/s; wind direction is from +155.0 to +92.0deg; drift angle is from +5.0 to 0.0deg. The variations of main flight variables for Aircraft B are also presented in Table 4.
The unsteady aerodynamic effects can be expected to be very significant under the circumstances of nearly rapid changes in angles of attack (α), flight altitude and Mach number at transonic flight. Since the highest α of Aircraft B reaches the values is above 6.6deg in transonic flight, compressibility effect is important. It should be noted that the turbulent vertical wind field was not measured or estimated in the QAR or FDR; but is included in the total α. The variation of wind speed for Aircraft B is the large, while drift angle is the small. The a z value of Aircraft A is not higher than Aircraft B, but the drop-off height and angles of attack variations are larger than those of Aircraft B.
3.2 Aircraft responses in operations
It can be observed that various flight variables change during severe clear-air turbulence encountered after analysing the aerodynamic and flight environments. Is it rule out that is the result of the aircraft operation to make the rapid changes of angles of attack (α), flight altitude and Mach number? It is of interest to examine the aircraft in operations in this section.
-
(1) Aircraft A
Figure 4 presents the variations of dynamic characteristics and control variables for Aircraft A at the altitude around 11,277m. The changes of α and sideslip angle ( $\beta $ ) including turbulence effects are shown in Fig. 4(a). The variation ranges of α are 6.2 to −5.9deg in turbulence encounter; the time history of $\beta $ is small, as indicated in Fig 4(a). The time history of $\theta $ and $\phi $ is shown in Fig. 4(b). The $\theta $ does not vary as much as α, but the highest value reaches 5.1deg The variations of $\phi $ is large with the variations −12 to 12deg during the turbulence encounter. The magnitudes of q and p are shown in Fig. 4(c); the variation of p is much larger than that of q.
The flight data recorder doesn’t show that the autopilot is deactivated during the turbulence encounter. The time history of $\delta_{\text{e}} $ and $\delta_{\text{a}} $ is shown in Fig. 4(d). When the altitude is dropping fast, $\delta_{\text{e}} $ becomes more negative, as shown in Fig. 4(d). On the other hand, around t = 6,855, when the aircraft is rising, $\delta_{\text{e}} $ becomes more positive. The variation of $\delta_{\text{e}} $ significantly affects the pitching moment coefficient and magnitude of q. It has no corresponding larger variations of elevator activity with respect to the variations of angles of attack. The values of $\phi $ are changed from negative to positive (t = 6,837–6,863s) and $\phi $ reaches a range from −12 to +12deg. While the $\phi $ is changing fast in the negative direction before t = 6,840s, $\delta_{\text{a}} $ is increased to positive, reaching +3deg at t = 6,845s, as shown in Fig. 4(d).
-
(2) Aircraft B
Figure 5 presents the variations of dynamic characteristics and control variables for Aircraft B at the altitude around 12,192m. The variation ranges of α are 6.6 to −3.4deg in turbulence encounter; the time history of $\beta $ is approach 0deg with small fluctuation, as indicated in Fig 5(a). The time history of $\theta $ and $\phi $ is shown in Fig. 5(b). The $\theta $ does not vary as much as α, but it has the magnitude with fluctuation. The variations of $\phi $ is large and the most magnitudes are in negative. The magnitudes of both data are fluctuation in abnormal appearance, it may cause due to the problems of measured sensors. The variables of p and q are unavailable in the FDR; the magnitudes of these two variables are estimated through the compatibility check. The magnitudes of both data are also fluctuation in abnormal appearance, as shown in Fig. 5(c). The flight data recorder doesn’t show that the autopilot is deactivated during the turbulence encounter. The time history of $\delta_{\text{e}} $ and $\delta_{\text{a}} $ is shown in Fig. 5(d).
The FDR of Aircraft A and B didn’t indicate that the autopilots were deactivated during the turbulence encounter. Based on the comprehensive analysis of the above two aircraft, the tentative conclusion is as follow:
-
(1) The magnitude variations of $\phi $ , p, and sideslip angle $\beta $ of Aircraft A are larger than those of Aircraft B.
-
(2) The magnitudes of both $\theta $ and $\phi $ of Aircraft B are fluctuation in abnormal appearance, these may be caused by the measured sensors.
-
(3) The results of aircraft operations causing drastic changes in angle of attack (α), flight altitude and Mach number have been ruled out.
3.3 Analysis of model predictions
In the present study, the accuracy of the established unsteady rolling moment model 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 coefficients and derivatives of pitching moment (C m) and rolling moment (C l) in the study of flight dynamic characteristics are calculated with these two models.
Figure 6 presents the C m and C l of the Aircraft A predicted by the unsteady pitching moment and rolling moment models, respectively. The predicted results by the final models have good match with the input data before the modeling. Once the C m and C l models are set up, one can calculate all necessary derivatives by central difference scheme to analyse the flight dynamic characteristics.
The final pitching- and rolling- moment models of Aircraft A consist of many fuzzy rules for each coefficient as described in Tables 5 and 6. In these tables, the numbers below each input variable represent the number of membership function. The total number of fuzzy cells (n) in each model is the product of each number, which is presented in column 3. The last column shows the final multiple correlation coefficients (R 2). The accuracy of the established pitching- and rolling-moment models through the fuzzy-logic algorithm can be judged by the (R 2).
Notes: Total number of fuzzy cells (n) = n 1 × n 2 × n3 × … n 10.
Notes: Total number of fuzzy cells (n) = n 1 × n 2 × n 3 × … n 11.
3.4 Influence of crosswind on loss of control
Aircraft A is subjected to crosswind before the plunging motion; the magnitudes of crosswind have the large variations in the time span of 6,830–6,860s, as shown in Fig. 7(a). The magnitude of drift angles is in the range of +5.0–0.0deg. The rolling motion in varying crosswind and the corresponding aileron deflections are presented in Fig. 7(b) and (c). The magnitudes of roll angles $\left( \phi \right)$ are changed from positive to negative (t = 6,835–6,855s) and $\phi $ reaches ±12deg in the time span of 6,838–6,850s. in Fig. 7(b), when the crosswind is abruptly decreased. While $\phi $ is changing fast in the time span of 6,830–6,880s i, the magnitudes of aileron angle ( $\delta_a $ ) are in large variations, as shown in Fig. 7(c). It is shown that the corresponding aileron angle will not be effective to control it. The roll control power derivatives should be positive in the effectiveness criteria of Table 3. Some values of $C_{l\delta a} $ are negative in plunging motion as shown in Fig. 7(d), and it implies uncontrollability because the corresponding aileron angle is ineffective. Although there are considerable control deflections of aileron angle ( $\delta_a $ ) around the time of fast variations of roll angles $\left( \phi \right)$ , and there is no control effectiveness in the period of t = 6,830–6,880s, it may be concluded that the rolling movement does not have well response to the aileron control due to the airflow detach with the control surface in a transient time.
The large variations of roll angles $\left( \phi \right)$ obviously induced by the influenced of varying crosswind. The corresponding aileron angle will not be effective to control the roll motion. It is of interest to examine the roll damping in oscillatory motion to evaluate the dynamic stability characteristics during the plunging motion.
Figure 8 presents the time history of roll damping in oscillatory motion for Aircraft A along the flight path to associate with $\dot \beta $ -derivatives. Note that in Fig. 8(c), the roll damping oscillatory derivatives are defined as:
The values of oscillatory derivatives are equivalent to the combinations of static damping and dynamic derivatives in above Equation (41). The use of oscillation derivative instead of static damping only is more consistent with the actual case in the analysis of dynamic stability characteristics. To be stable, (C lp ) osc < 0. Physically, if it is unstable, the motion will be divergent in oscillatory motions. The values in the period of plunging motion have some differences between oscillatory and static damping derivatives in Fig. 8(a) and (c) for Aircraft A due to the effects of the dynamic derivatives (i.e. $\dot \beta $ -derivatives).
Figure 8(a) and (b) show roll damping and dynamic derivatives of stability, respectively. To be stable, the roll damping (C lp < 0) and dynamic derivatives of stability ( ${C_{l\dot \beta }}$ < 0). The magnitudes of C lp are positive with insufficient in roll damping and ${C_{l\dot \beta }}$ are positive with unstable conditions in the period of t = 6,835–6,860s. The effects of $\;\dot \beta $ -derivative on (C lp ) osc are to cause the lateral-stability more unstable in Fig. 8(c) through the magnitude comparisons of (C lp ) osc with C lp . It implies that the effects of $\;\dot \beta $ -derivativecan cause the lateral-stability more unstable.
3.5 Sources of significant increase in angles of attack
Figure 9 presents the time history of flight variables and with pitch-control-power derivative for Aircraft A during cruise phase. The flight variables in Fig. 9(a) are vertical acceleration (a z ), angle-of-attack (α), pitch angle ( $\theta $ ), elevator angle ( $\delta_e $ ). The variations of a z and α show the highest magnitudes being 1.67g and 6.2deg at t = 6,863s and the lowest being close to 0g and 5.9deg around t = 6,856s, respectively. Figure 9(b) shows that variations of α and $\theta $ are approximately in phase with a z during turbulence encounter. The magnitude of α should be approximately close to $\theta $ , if there is no airflow disturbance during cruise flight. The magnitudes of α are obviously larger than those of $\theta $ during turbulence encounter in the time span between t = 6,830–6,865s.
At transonic speeds, the variation range of angle-of-attack is +6.2 to −5.9deg as shown in Fig. 9(b) is high for transport aircraft in cruise. Numerical calculation of a typical supercritical aerofoil shows strong shock and possible flow separation at 6.5deg in angle-of-attack. Yet on the subject aircraft, the largest angle-of-attack reached is 6.2deg. Is this large angle produced by the elevators? The elevator angles ( $\delta_e $ ) and the corresponding control power derivatives ( ${C_{{m_{\delta e}}}}$ < 0) are presented in Fig. 9(c) and (d), respectively. Under normal conditions, the pitch control power derivatives should be negative ( ${C_{{m_{\delta e}}}}$ < 0), meaning the pitching moment being negative if the elevator being positive (i.e. trailing-edge down). From Fig. 9(d), it is seen that fast angle-of-attack variation during the period t = 6,830–6,880s, the longitudinal control of $C_{m\delta e} $ is ineffective. It may be concluded that the pitching movement does not have a good response to the elevator control due to the airflow detach with the control surface in a transient time.
The most parts of elevator activities are nose-up, because they are negative in elevator angles as shown in Fig. 9(c); except the activity of elevators at t = 6,855–6,865s. is nose-down. In general, the magnitudes of angles of attack and pitch angles should be increased with the equivalent values in cruise flight as the response respect to nose-up elevator activities. However, the most parts of angles of attack are larger than pitch angles as shown in Fig. 9(b) and some parts of control power derivatives are positive as shown in Fig. 9(d) in this period. The variation of angle-of-attack is essential, but it has less corresponding variation of elevator activity in this period. Therefore, some increments in angle-of-attack are not produced by elevators.
Equation (1) indicates that mergers of α m and α w are the measurement of the alpha sensor due to the aircraft motion and wind speed influence (including vertical wind speed), respectively. The magnitudes of α are rapidly changed in phase with a z , therefore, the mutation of α should be mainly influenced by vertical wind speed. Are certain variations of α produced by the aircraft motion? The movement mechanism will be further investigations below.
The pitching moment (PM ic) represents pitch rate coupling with yaw rate in Equation (20). PM ic is very important in fighter design and not limited to high-speed conditions. If both “p” and “r” are of the same sign, and I zz > I xx, (usually being the case), PM ic is positive.
Figure 10(a), (b) and (c) presents the variations of roll rate, yaw rate and pitching moment due to inertial coupling (PM ic). PM ic has a large variation during t = 6,830–6,860s. due to the large variations of roll and yaw rates, especially for roll rate. The rapid variations of angles of attack should be contributed by PM ic during this period. A tentative conclusion is that if the variations of angle-of-attack and pitch angle occur at the same time, the increase or decrease must be caused by a pitching moment, or inertial coupling in this case. This is because the pitch angle can only be changed by a pitching moment. On the other hand, if the pitch angle change is not in phase with that of the angle-of-attack, the increase in angle-of-attack should be caused by upflow in the turbulent air, which is usually assumed in random motion [Reference Lee and Lan23, Reference Sheu and Lan24], resulting in flight altitude loss.
4.0 Concluding remarks
The objective of the present paper was to present the movement mechanisms of transport aircraft response to severe clear-air turbulence to obtain the loss of control prevention for pilot training in IATA – Loss of Control In-flight (LOC-I) program. The comparative analyses of flight environment and aircraft response to severe clear-air turbulence for two four-jet aircraft were studied. The one with larger dropped-off altitude during the plunging motion was chosen to construct the movement mechanism. The nonlinear unsteady aerodynamic model of the chosen transport aircraft was established through flight data mining and the fuzzy-logic modeling techniques based on post-flight data. The influences of varying vertical wind and crosswind on loss of control were presented. The source of significant angles of attack induced by the varying crosswind during the sudden plunging motion were proved through the analysis of oscillatory derivatives based on model predicted data. The numerical results and discussions of application to the movement mechanisms were concluded as follows:
-
(1) The crosswind before the turbulence encounter would easily induce rolling motion and then initiate the sudden plunging motion during the turbulence encounter. The magnitude of varying crosswind encounter on the flight route could be judged by the drift angle.
-
(2) The roll rate would increase the oscillatory rolling motion during plunging motion, if the rolling damping was insufficient. The drop-off altitude would be enlarged by the oscillatory rolling motion during the sudden plunging motion.
-
(3) The variation of angle-of-attack is essential, but it has less corresponding variation of elevator activity in severe clear-air turbulence. The increments of α should be mainly influenced by vertical wind speed; some increments in angle-of-attack were caused by the pitching moment due to inertial coupling PMic after the movement mechanism being further investigated.
The analytical results could provide the mitigation concepts and promote the understanding of dynamic responses of the transport aircraft in severe clear-air turbulence. To formulate preventive actions, the situation awareness of drift angle variations for the operational pilot before and during severe atmospheric turbulence encounter will be the further study task in the future. The present concept is an innovation to improve aviation safety and operational efficiency for the airlines. It is expected to provide a valuable lecture for international training courses for the IATA – Loss of Control In-flight (LOC-I) program after this paper is published.
Data and information availability statement
The flight data of FDR from the specific aviation accident to support the present study is restricted to deliver to the third party. The flight number, name of airlines and aircraft manufacture are not publicly available without permission.
Acknowledgements
This research project is sponsored by the grant of 2020M671385 from China Postdoctoral Science Foundation and BK20200175 of Jiangsu basic research program (NSFC) – Youth Fund Project.