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.

Figure 1. Flowchart for development of a non-linear unsteady aerodynamic database.

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.

Table 1. Raw data of wind filed

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.

(1) \begin{equation}{\rm{Angle\ of\ attack}}\ \alpha = {\alpha _{\rm{m}}} + {\alpha _{\rm{w}}}\end{equation}

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]:

(2) \begin{equation}\dot \phi = p + q\sin \phi \tan \theta + r\cos \phi \tan \theta \end{equation}

(3) \begin{equation}\dot \theta = q\cos \phi - r\sin \phi \end{equation}

(4) \begin{equation}\dot \psi = \left( {q\sin \phi + r\cos \phi } \right)\sec \theta \end{equation}

(5) \begin{equation}\dot V = \left({a_x} - g\sin \theta \right)\cos \alpha \cos \beta + \left( {{a_y} + g\sin \phi \cos \theta } \right)\sin \beta + \left( {{a_z} + g\cos \phi \cos \theta } \right)\sin \alpha \sin \beta \end{equation}

(6) \begin{equation}\dot \alpha = \frac{{\left[ {\left( {{a_z} + g\cos \theta \cos \phi } \right)\!\cos \alpha - \left( {{a_x} - g\sin \theta } \right)\sin \alpha } \right]}}{{(V\cos \beta )}} + q - \tan \beta \!\left( {p\cos \alpha + r\sin \alpha } \right)\end{equation}

(7) \begin{align}\dot \beta & = \frac{{\cos \beta \!\left( {{a_y} + g\cos \theta \sin \phi } \right)}}{V} + p\sin \alpha - r\cos \alpha\nonumber\\& \quad - \sin \beta \!\left[\left({a_z} + g\cos \theta \cos \phi \right)\sin \alpha - ({a_x} - g\sin \theta )\cos \alpha \right]/V\end{align}

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:

(8) \begin{equation}\dot{\mathop Z\limits^{\rightharpoonup}} = {\mathop f\limits^{\rightharpoonup}} \!\left( {\mathop x\limits^{\rightharpoonup}} \right) = {\mathop f\limits^{\rightharpoonup}} \left( {{{\mathop x\limits^{\rightharpoonup}}_m} - \Delta {\mathop x\limits^{\rightharpoonup}} } \right)\end{equation}

(9) \begin{equation}{\mathop Z\limits^{\rightharpoonup}} = {\left( {V,\bar \alpha ,\bar \beta ,\theta ,\phi ,\psi } \right)^T}\end{equation}

(10) \begin{equation}{{\mathop x\limits^{\rightharpoonup}} _m} = {\left( {{{\bar a}_x},{{\bar a}_y},{{\bar a}_z},\bar p,\bar q,\bar r,V,\bar \alpha ,\bar \beta ,\theta ,\phi ,\psi } \right)^T}\end{equation}

(11) \begin{equation}\Delta {\mathop x\limits^{\rightharpoonup}} = {\left( {{b_{_x}},{b_{{a_y}}},{b_{{a_z}}},{b_p},{b_q},{b_r},{b_V},{b_\alpha },{b_\beta },{b_\theta },{b_\phi },{b_\psi }} \right)^T}\end{equation}

where overbar “–” stands for the mean value, and the subscript “m” indicates the measured or recorded values.

The cost function is defined as:

(12) \begin{equation}J = \frac{1}{2}{\left( {\dot{\mathop z\limits^{\rightharpoonup}} - {\mathop f\limits^{\rightharpoonup}} } \right)^T}Q\!\left( {\dot{\mathop z\limits^{\rightharpoonup}} - {\mathop f\limits^{\rightharpoonup}} } \right)\end{equation}

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.

Table 2. Main geometry and moment of inertia data for Aircraft A and B

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. (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. (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 (RM_ic) due to inertia coupling as follows:

(19) \begin{equation}R{M_{ic}} = qr\!\left( {{I_{zz}}-{I_{yy}}} \right)\end{equation}

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:

(20) \begin{equation}P{M_{ic}} = rp\!\left( {{I_{zz}}-{I_{xx}}} \right)\end{equation}

(21) \begin{equation}Y{M_{ic}} = pq\!\left( {{I_{yy}}-{I_{xx}}} \right)\end{equation}

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:

(22) \begin{equation}T = f\!\left( {h,W,M,\textrm{CAS},\textrm{EPR},{{\dot m}_f}} \right)\end{equation}

For General Electrical (GE) or CFM International (CFM) turbofan engines, the rpm of the low-pressure compressor (N ₁ ) is used to set the level of thrust, so that the thrust model is set up as:

(23) \begin{equation}T = f\!\left( {h,W,M,\textrm{CAS},{N_1},{{\dot m}_f}} \right)\end{equation}

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.

Figure 2. Flowchart for platform development of root cause summary in aviation accident.

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]:

(24) \begin{equation}{C_x},{C_z},{C_m} = f\!\left( {\alpha ,\dot \alpha ,q,{k_1},\beta ,{\delta _e},M,p,{\delta _s},\bar q} \right)\end{equation}

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 ₁ ), 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],

(25) \begin{equation}{C_y},{C_l},{C_n} = f\!\left( {\alpha,\beta ,\phi ,p,r,{k_2},{\delta _a},{\delta _r},M,\dot \alpha ,\dot \beta } \right)\end{equation}

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 ₂ ), 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):

(26) \begin{equation}f(x + \Delta x) = f(x) + f{\rm{'}}(x)\Delta x + \frac{{{f^{(2)}}(x)}}{2}\Delta {x^2} + \ldots \end{equation}

Taylor expansion to find f(x−Δx):

\begin{equation*}f(x - \Delta x) = f(x) + f{\rm{'}}(x)({-}\Delta x) + \frac{{{f^{(2)}}(x)}}{2}{({-}\Delta x)^2} + \ldots \end{equation*}

Then,

(27) \begin{equation}f(x - \Delta x) = f(x) - f{\rm{'}}(x)\Delta x + \frac{{{f^{(2)}}(x)}}{2}\Delta {x^2} + \ldots \end{equation}

Subtracting (26) from (28), one gets

(28) \begin{equation}f(x + \Delta x) - f(x - \Delta x) = 2f'(x)\Delta x\end{equation}

Therefore,

(29) \begin{equation}f'\left( x \right) = {\frac{f\!\left( {x + \Delta x} \right) - f\!\left( {x - \Delta x} \right)}{2\Delta x}}.\end{equation}

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:

(30) \begin{equation}{C_m}_a = [{C_m}(\alpha + \Delta \alpha ,{\rm{ }} \mbox{- - -} ) - {C_m}(\alpha - \Delta \alpha ,\mbox{- - -} )]/2\Delta \alpha \end{equation}

Δα = 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:

(31) \begin{equation}{C_{l{\rm{ }}p}} = [{C_l}( \mbox{- - -} ,p\Delta + p, \mbox{- - -} ) - {C_l}( \mbox{- - -} ,p - \Delta p, \mbox{- - -} )]/2\Delta p\end{equation}

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.

Table 3. The effectiveness criteria for derivative values of stability, damping and controllability [Reference Roskam16]

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]:

(32) \begin{equation}{P^i} = {y_i}\!\left( {{x_1},{x_2}, \cdots\! ,{x_r}, \cdots {x_k}} \right) = P_0^i + P_1^i{x_1} + \cdots + P_r^i{x_r} + \cdots P_k^i{x_k}\end{equation}

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

(33) \begin{equation}{x_{r,norm}} = \frac{{{x_r} - {x_{r,min}}}}{{{x_{r,max}} - {x_{r,min}}}}\;,r = 1,2, \cdots\!,k\end{equation}

Figure 3. Flowchart of refining process of fuzzy-logic modeling.

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:

(34) \begin{equation}\hat y = \frac{{\mathop \sum \limits_{i = 1}^n op\!\left[ {{A^i}\!\left( {{x_{1,j}}} \right), \cdots\!,{A^i}\!\left( {{x_{r,j}}} \right), \cdots\!,{A^i}\!\left( {{x_{k,j}}} \right)} \right]{P^i}}}{{\mathop \sum \limits_{i = 1}^n op\!\left[ {{A^i}\!\left( {{x_{1,j}}} \right), \cdots\!,{A^i}\!\left( {{x_{r,j}}} \right), \cdots\!,{A^i}\!\left( {{x_{k,j}}} \right)} \right]}}\;j = 1,2, \cdots\!,m\end{equation}

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]:

(35) \begin{equation}SSE = \mathop \sum \limits_{j = 1}^m\! {\left( {{{\hat y}_j} - {y_j}} \right)^2}\end{equation}

On the other hand, the structure of fuzzy cells is optimised by maximising the multiple correlation coefficients (R ² ):

(36) \begin{equation}{R^2} = 1 - \frac{{\left\{ {\mathop \sum \limits_{j = 1}^m {{\left( {{{\hat y}_j} - {y_j}} \right)}^2}} \right\}}}{{\left\{ {\mathop \sum \limits_{j = 1}^m {{\left( {\bar y - {y_j}} \right)}^2}} \right\}}}\end{equation}

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 ⁱ –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:

(37) \begin{equation}P_{r,t + 1}^i = P_{r,t}^i - {\alpha _r}\frac{{\partial \!\left( {SSE} \right)}}{{\partial P_r^i}}\end{equation}

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

(38) \begin{equation}\frac{{\partial \!\left( {SSE} \right)}}{{\partial P_r^i}} = 2\mathop \sum \limits_{j = 1}^m \!\left( {{{\hat y}_j} - {y_j}} \right)\frac{{\partial {{\hat y}_j}\left( {{x_{1,j}}, \cdots\!,{x_{k,j}},P_r^1, \cdots\!,P_k^n} \right)}}{{\partial P_r^i}}\end{equation}

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,

(39) \begin{equation}P_{0,t + 1}^i = P_{0,t}^i - 2{\alpha _0}\!\left( {{{\hat y}_j} - {y_j}} \right)\frac{{op\!\left[ {A_1^i\!\left( {{x_{1,j}}} \right), \cdots\!,A_k^i\!\left( {{x_{k,j}}} \right)} \right]}}{{\mathop \sum \limits_{s = 1}^n op\!\left[ {A_1^s\!\left( {{x_{1,j}}} \right), \cdots\!,A_k^s\!\left( {{x_{k,j}}} \right)} \right]}}\end{equation}

and for r = 1, …, k,

(40) \begin{equation}P_{r,t + 1}^i = P_{r,t}^i - 2{\alpha _r}\!\left( {{{\hat y}_j} - {y_j}} \right)\frac{{op\!\left[ {A_1^i\!\left( {{x_{1,j}}} \right), \cdots\!,A_k^i\!\left( {{x_{k,j}}} \right)} \right]{x_{r,j}}}}{{\mathop \sum \limits_{s = 1}^n op\!\left[ {A_1^s\!\left( {{x_{1,j}}} \right), \cdots\!,A_k^s\!\left( {{x_{k,j}}} \right)} \right]}}\end{equation}

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. (1) Cost = SSE _t = minimised

2. (2) $RER = \dfrac{{SS{E_t} - SS{E_{t - 1}}}}{{SS{E_t}}} \lt {\varepsilon _2}$

3. (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.

Table 4. The variations of main flight variables in aerodynamic and flight environments

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. (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.

Figure 4. Variations of dynamic characteristics and control variables for Aircraft A.

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).

1. (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).

Figure 5. Variations of dynamic characteristics and control variables for Aircraft B.

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. (1) The magnitude variations of $\phi $ , p, and sideslip angle $\beta $ of Aircraft A are larger than those of Aircraft B.

2. (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. (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 ²). 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 ²). The accuracy of the established pitching- and rolling-moment models through the fuzzy-logic algorithm can be judged by the (R ²).

Table 5. Final models of pitching-moment coefficients of Aircraft A

Notes: Total number of fuzzy cells (n) = n ₁ × n ₂ × n₃ × … n ₁₀.

Table 6. Final models of rolling-moment coefficients of Aircraft A

Notes: Total number of fuzzy cells (n) = n ₁ × n ₂ × n ₃ × … n ₁₁.

Figure 6. The comparisons of the predicted results and input data before the modeling for Aircraft A.

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.

Figure 7. Influence of varying crosswind corresponding to rolling control for Aircraft A.

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:

(41) \begin{equation}{\left( {{C_{lp}}} \right)_{osc}} = {C_{lp}} + {C_{l\dot \beta }}\sin \alpha \end{equation}

Figure 8. The time history of roll damping and oscillatory derivatives for Aircraft A.

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.

Figure 9. Time history of flight variables and with pitch-control-power derivative for Aircraft A.

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.

Figure 10. Roll rate, yaw rate and pitching moment due to inertial coupling.