# Fuzzy Logic-Based Aerodynamic Modeling with Continuous Differentiability

Authors: Ray C. Chang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper presents a modeling method based on a fuzzy-logic algorithm to establish aerodynamic models by using the datasets from flight data recorder (FDR). The fuzzy-logic aerodynamic models are utilized to estimate more accurately the nonlinear unsteady aerodynamics for a transport aircraft, including the effects of atmospheric turbulence. The main objective in this paper is to present the model development and the resulting models with continuous differentiability. The uncertainty and correlation of the data points are estimated and improved by monitoring a multivariable correlation coefficient in the modeling process. The latter is increased by applying a least square method to a set of data points to train a set of modeling coefficients. A commercial transport aircraft encountered severe atmospheric turbulence twice at transonic flight in descending phase is the study case in the present paper. The robustness and nonlinear interpolation capability of the fuzzy-logic algorithm are demonstrated in predicting the degradation in performance and stability characteristics of this transport in severe atmospheric turbulence with sudden plunging motion.


---

## Body

## 1. Introduction

During the operation of aircraft, it is very difficult to avoid the influence of adverse weather conditions. The effects of adverse atmospheric conditions usually are related to the dynamic aerodynamic effects that result from instantaneous changes of aircraft flight attitude. This situation, which cannot be easily flight-tested in the aircraft certification process, is the so-called abnormal flight condition [1]. Such abnormal flight condition is very likely to cause the deployed aircraft out of control. In recent years, the transport aircraft response to the hazardous weathers due to atmospheric disturbance has been of great concern. The windshear, down burst, and atmospheric turbulence of hazardous weathers often lead to a high level of dynamic aerodynamic effects. These hazardous weathers usually have the characteristics of strong crosswind or vertical gust. It is not possible with the current technology to simulate the atmospheric turbulence effects by conventional flight simulators or wind tunnel on the aerodynamics, stability, and flying quality of aircraft [2]. In addition, the aerodynamic derivative estimation from flight test data is always based on the concept of a small disturbance theory around an equilibrium state without having the influences of dynamic aerodynamic effects.To provide the mitigation concepts and promote the understanding of aerodynamic responses of the transport aircraft under adverse weather conditions, a research project to evaluate the effects of weather hazards on performance and stability characteristics for transport aircraft is undertaken. To effectively analyze the performance degradation and variations in static and dynamic stability of transport aircraft encountering hazardous weather, the nonlinear and dynamic (i.e., time-dependent) aerodynamic models based on flight data would be needed.Traditional approaches of flight data analysis, such as the maximum likelihood method (MMLE) [3], the least-square or the regression method [4], apply the identification methods to data over an interval of time or other state variables to identify models with a quasisteady assumption over the interval. However, these conventional methods have not been shown to be capable of handling nonlinear and unsteady (i.e., dynamic) aerodynamic environments possibly exhibited in abnormal flight condition, in particular in hazardous weather, such as in-flight icing, windshear, and atmospheric turbulence.In setting up the nonlinear and unsteady aerodynamic models to predict the required continuous derivatives of the function in aerodynamics for the analyses of stability and controllability, the traditional methods of using the tabulated data and the performance figures are very difficult to correlate complex functional relations among numerous parameters.In earlier development of the fuzzy-logic algorithm, Zadeh [5] used the fuzzy sets to simulate physical parameters with membership functions. The disadvantage in this approach is that the predicted curves in the functional approximation tend to be piecewise continuous. Therefore, it was not suitable when derivatives of the function were needed in some applications of aviation technology. In 1985, Takagi and Sugeno [6] used the internal functions, instead of the fuzzy sets, in developing the fuzzy-logic algorithm.It is well known that fuzzy models proposed by Takagi and Sugeno can be effectively represent unsteady and nonlinear system. Many important results in the applied areas, such as sensor fault detection [7], stabilization systems with time-varying delay [8], discrete-time state-delay systems [9], nonlinear itô stochastic systems with time-varying delay [10], a flexible air-breathing hypersonic vehicle [11], and fuzzy stochastic systems with time-varying delay [12] have great efforts in modern science and engineering, especially for some systems without having sufficient information as input databases for modeling [13, 14].In this approach, the main features are the internal functions, membership functions, and the output cells. The output curves in the prediction are smooth. Tan and Xie [15, 16] applied the theory to simulate microelectronic processes with good accuracy. The present paper is based on the modeling technique suggested by Tan and Xie. This technique was first applied to aviation technology in 1997 and later to flight data [17–24]. A new study to examine the effects of severe atmospheric turbulence on aircraft performance and stability in descending flight is presented in this paper.Atmospheric turbulence has long been a difficult issue for the aviation community because it is not only the leading cause of serious injury in nonfatal accidents, but also most likely to cause the aircraft loss of control. Both the Federal Aviation Administration and National Aeronautics and Space Administration have research programs to deal with this problem. One of the objectives in NASA’s Aviation Safety Program is to study the turbulence hazards through research, flight experiments, and data analysis. The hazard levels are quantified in terms of the peak values of the root mean square (RMS) of the normal loads with moving 5-sec intervals [25, 26]. The estimation was based on the assumption of von Karman’s power spectral density [27]. When von Karman's power spectral density is assumed, the Gaussian distribution is utilized to exhibit continuous turbulence. However, when localized severe turbulence or gust is included, the results will fall way off the correlation curve for the peak RMS normal loads to peak loads. Experience indicates that one main type of motion to cause flight injuries in atmospheric turbulence is the sudden plunging motion with the abrupt change in altitude.A commercial transport aircraft encountered severe atmospheric turbulence twice in the descending phase. As a result, several passengers and cabin crews sustained injuries, because of which this event was classified as the aviation accident. Since one main type of motion to cause flight injuries is the sudden plunging motion with the abrupt change in altitude at transonic flight, so the main objective of this paper is to present the model development and the capability of resulting models to examine the static and dynamic stability characteristics in severe atmospheric turbulence with sudden plunging motion.


## 2. Theoretical Development

Since the nonlinear unsteady 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 system model and the method to identify parameters of 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. In order to obtain all the corresponding output, the present paper uses internal functions instead of fuzzy sets [5, 28, 29] to generate the output of the model.The general idea of the present FLM technique is to set up the relations between its input and output variables of the whole system. The internal functions, membership functions, and outputs are three basic elements for the FLM approach. Two main tasks are involved in the FLM process. One is the identification of the coefficients of the internal functions, which is called parameter identification. The other one is structure identification to identify the optimal structure of fuzzy cells of the model. Details of the FLM technique are described in the followings.


### 2.1. Internal Functions

The fuzzy-logic model uses many internal functions to cover the defined ranges of the influencing parameters. Although the form of internal functions is very simple, it can represent a highly nonlinear relationship between the input and output for the whole system. These internal functions are assumed to be linear functions of input variables [15, 16] as follows:
(1)Pi=yi(x1,x2,…,xr,…,xk)=p0i+p1ix1+⋯+prixr+⋯+pkixk,
where pri, r=0,1,2,…,k, are the coefficients of internal functions yi, and k is number of input variables.The recorded data in QAR, such as flight altitude (h), calibrated airspeed (CAS), and angle of attack (α), is chosen as the input variables for a specific fuzzy model. In the present paper, yi is denoted as an estimated aerodynamic coefficient of force or moment, and xr are the variables of the input data. The numbers of the internal functions (i.e., cell’s numbers) are quantified by the membership functions.


### 2.2. 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 withA(xr) as its membership grade. One membership function from each variable constitutes a fuzzy cell. For the ith cell, the corresponding membership grades are represented by Ari(xr), 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 calculated. For a given system with input variables x1,x2,…,xr,…,xk, the recorded values of each input variables are normalized by using (xr-xr,min)/(xr,max-xr,min) to transform them into the ranges of [0,1]. The membership grading also ranges from 0 to 1; “0” meaning no effect from the corresponding internal function, and “1” meaning a full effect. Generally, overlapped straight lines, triangles, or trapezoids are frequently the shapes used to represent the functions. Because overlapped triangular membership function is simple and involves less computing time, it represents the grades of membership functions in the present FLM technique.The total number of fuzzy cells isn=N1×N2×⋯×Nr×⋯×Nk. For a variable xr, the number of membership function is Nr. Each fuzzy cell is in a different combination from others formed by taking one membership function from each input variable.In the present application, triangular membership functions are used throughout. LetN be the number of membership functions and j the index for thejth membership functions. Then the membership grades can be described as follows:(1)
N=2:
(2)A(xr)=xr,j=1,A(xr)=1-xr,j=2,(2)
N≥3:forj=3 to N-m, where m is equal to the greater number of 0 and integer of (N-2)/2:
(3)A(xr)=xrdu,0≤xr≤du,A(xr)=(1-xr)(1-du),du≤xr≤1,wheredu=Δx1*(j-2), and Δx1=1.0/(N-m-1).Forj≥N-m(4)A(xr)=(dd-xr)dd,0≤xr≤dd,A(xr)=(dd-xr)(dd-1),dd≤xr≤1,wheredd=Δx2*(j-N+m), and Δx2=1.0/(m+1).The membership functions of triangular shapes are illustrated in Figure1. In Figure 1, although the membership functions are continuous functions, there are discontinuities in slopes at some points. However, differentiation of membership functions is not performed in estimating derivatives. In the present application, aerodynamic derivatives are all estimated with a central difference scheme, which will be presented later.Triangular membership functions.
(a)
(b)
(c)
(d)


### 2.3. Fuzzy Rule Inference

A fuzzy cell is formed by taking one membership function from each variable. The total number of cells is the number of possible combinations by taking one membership function from each input variable. For every cell, it has a fuzzy rule to guide the input and output relations. The rule of theith cell [15, 16] is stated as
(5)ifx1isA1i(x1),andx2isA2i(x2),and…,andxkisAki(xk),thenthecelloutputisPi=p0i+p1ix1+⋯+prixr+⋯+pkixk,
where i=1,2,…,n is the index of the cells, n is the total number of cells of the model; Pi(x1,x2,…,xr,…,xk) is the internal function with parameters p0i,p1i,…,pri,…,pki to be determined, and Aki(xk) denotes the membership function for xk. Each function covers a certain range of input variables.


### 2.4. Defuzzification

In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function, (5). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm (i.e., fuzzy rule inference). Because of this weighting among many factors over large ranges of possibilities, the word “fuzzy” is derived to describe the method. However, its prediction after defuzzification is in a deterministic value (a ‘‘crisp’’ one). The output estimated by the fuzzy-logic algorithm corresponding to the jth input (x1,j,x2,j,…,xr,j,…,xk,j) is as follows:
(6)y^j=∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]pi∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)].
In (6) product [Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)] is the weighted value of the ith 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 “product” stands product operator of its elements in this paper.


### 2.5. Parameter Identification

Given a set of membership functions for each input variable, the unknown coefficients of the internal functions are determined by using the gradient-descent method. The gradient descent is also known as steepest descent. The accuracy of the established aerodynamic model through the fuzzy-logic algorithm is estimated by the sum of squared errors (SSEs) and the multiple correlation coefficients(R2):
(7)SSE=∑j=1m(y^j-yj)2,(8)R2=1-{∑j=1m(y^j-yj)2}{∑j=1m(y--yj)2}.
In (7) and (8), where y^j, the output of the fuzzy-logic model at point j, is estimated by (6); yj is the data point used for the model training at point j; y- is the mean of the sample data, and m is the total number of data points. The model training is to determine the unknown coefficients of the internal functions, pri, by maximizing the value of R2. These coefficients are determined by the following iterative formula to minimize the sum of squared error (see (7)):
(9)pr,t+1i=pr,ti-αr∂(SSE)∂pri,
where αr is the convergence factor or the step size in the gradient method; subscript index t denotes the iteration sequence.After simplification, (9) becomes the following.Forr=0,(10a)p0,t+1i=p0,ti-2α0(y^j-yj)=×product[A1i(x1),…,Ari(xr,j),…,Aki(xk,j)]∑s=1nproduct[A1s(x1,j),…,Ars(xr,j),…,Aks(xk,j)],
and for r=1,2,…,k,
(10b)pr,t+1i=pr,ti-2αr(y^j-yj)=×product[A1i(x1),…,Ari(xr,j),…,Aki(xk,j)]xr,j∑s=1nproduct[A1s(x1,j),…,Ars(xr,j),…,Aks(xk,j)].The present FLM technique was utilized to conduct NASA research projects [17, 18]. The unsteady aerodynamic models for F-16XL were based on dynamic wind-tunnel data in [17], and those were from flight test data in [18]. The longitudinal and lateral-directional derivatives would be useful in aircraft stability analysis and control system design. Although the data source or application was not the same, the iteration during the search sequence stops of modeling should be the same with a principle for judgment.The iteration during the search sequence stops when one of the following three criteria is satisfied:(11)(1)Cost=SSEt<ε1,(2)RER=SSEt-SSEt-1SSEt<ε2,(3)t=tmax.
In the previos criteria, Cost=SSEI is the sum of squared errors (SSE) in current iteration to be denoted by “Cost” and RER=(cost_current-cost_previous)/cost_current to be denoted by “RER” (i.e., the relative error) for simplicity in descriptions; ε1 and ε2 are the required precision criteria; and tmax is a specified maximum iteration numbers.Given membership functions and the training data, this parameter identification procedure can be applied to establish a fuzzy-logic model. Although the procedure is well understood, to obtain a good model with the appropriate membership functions is quite challenging and requires numerical experimentation.


### 2.6. Model Structure Identification

In the fuzzy-logic model, the model structure is indicated by the number of membership functions for each variable. For a fuzzy-logic model with multiple variables, the structure is the combination of the numbers and forms of the membership functions assigned to all input variables. Since the sequence defines the one-to-one relationship between the numbers and the forms for each variable, the structure can be uniquely described by numbers of the membership functions.The model structure is determined by maximizing the correlation coefficient, (8). A search forward algorithm has been employed for the identification. At each search stage, there may be many fuzzy-logic models with different structure combinations. The search stage numbers are denoted by Ns. Out of all the possible intermediate fuzzy-logic models at each search stage, for an efficient search, only some structures are developed and evaluated. Two selection criteria, to be given below, are used to choose these structures. With the incremental sequence and the selection criteria, the search forward algorithm is summarized as follows.(1)
Specify the input variablesxr, r=1,2,…,k and the output variable y.(2)
Assume an initial structure also called parent structure as(N10,N20,…,Nr0,…,Nk0).(3)
Begin at the search stage numberNs=1, form all possible structures starting from the parent structure by adding one more membership function a time only to one input variable. Those all possible structures are called child structures as (N10+1,N20,…,Nr0,…,Nk0), (N10,N20+1,…,Nr0,…,Nk0),…,(N10,N20,…,Nr0,…,Nk0+1). Perform the identification of internal coefficients in (1) for each child structure and then calculate the R2 by using (8).(4)
Select the top 5 child structures among all calculated values ofR2 as new parent structures for next search step Ns=Ns+1.(5)
Go back to step(2) starting from the new parent structures and repeat the same procedures in steps (2) and (3) until the best structure is identified.(6)
Pick out the maximum value ofR2 among the child structures in each searching stage as Rmax2. The structure with the largest Rmax2 corresponding to all picking values is the optimal structure within a sensible Ns.The previous process is illustrated in Figure2. The structure identification is tied up with parameter identification. In the structure identification, parameter identification to determine the p-parameters according to (9) is also needed. The value of Nr is gradually changed in each iteration until the composite of Nr (i.e., optimal model structure) being fund, while the searching stage reaches Rmax2, for the model structure identification; and the composite of Nr is fixed, the number of iteration is unfixed until both the values of R2 and RER reach the requirements in the final the parameter identification.Figure 2
Identification process for the best structure.


## 2.1. Internal Functions

The fuzzy-logic model uses many internal functions to cover the defined ranges of the influencing parameters. Although the form of internal functions is very simple, it can represent a highly nonlinear relationship between the input and output for the whole system. These internal functions are assumed to be linear functions of input variables [15, 16] as follows:
(1)Pi=yi(x1,x2,…,xr,…,xk)=p0i+p1ix1+⋯+prixr+⋯+pkixk,
where pri, r=0,1,2,…,k, are the coefficients of internal functions yi, and k is number of input variables.The recorded data in QAR, such as flight altitude (h), calibrated airspeed (CAS), and angle of attack (α), is chosen as the input variables for a specific fuzzy model. In the present paper, yi is denoted as an estimated aerodynamic coefficient of force or moment, and xr are the variables of the input data. The numbers of the internal functions (i.e., cell’s numbers) are quantified by the membership functions.


## 2.2. 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 withA(xr) as its membership grade. One membership function from each variable constitutes a fuzzy cell. For the ith cell, the corresponding membership grades are represented by Ari(xr), 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 calculated. For a given system with input variables x1,x2,…,xr,…,xk, the recorded values of each input variables are normalized by using (xr-xr,min)/(xr,max-xr,min) to transform them into the ranges of [0,1]. The membership grading also ranges from 0 to 1; “0” meaning no effect from the corresponding internal function, and “1” meaning a full effect. Generally, overlapped straight lines, triangles, or trapezoids are frequently the shapes used to represent the functions. Because overlapped triangular membership function is simple and involves less computing time, it represents the grades of membership functions in the present FLM technique.The total number of fuzzy cells isn=N1×N2×⋯×Nr×⋯×Nk. For a variable xr, the number of membership function is Nr. Each fuzzy cell is in a different combination from others formed by taking one membership function from each input variable.In the present application, triangular membership functions are used throughout. LetN be the number of membership functions and j the index for thejth membership functions. Then the membership grades can be described as follows:(1)
N=2:
(2)A(xr)=xr,j=1,A(xr)=1-xr,j=2,(2)
N≥3:forj=3 to N-m, where m is equal to the greater number of 0 and integer of (N-2)/2:
(3)A(xr)=xrdu,0≤xr≤du,A(xr)=(1-xr)(1-du),du≤xr≤1,wheredu=Δx1*(j-2), and Δx1=1.0/(N-m-1).Forj≥N-m(4)A(xr)=(dd-xr)dd,0≤xr≤dd,A(xr)=(dd-xr)(dd-1),dd≤xr≤1,wheredd=Δx2*(j-N+m), and Δx2=1.0/(m+1).The membership functions of triangular shapes are illustrated in Figure1. In Figure 1, although the membership functions are continuous functions, there are discontinuities in slopes at some points. However, differentiation of membership functions is not performed in estimating derivatives. In the present application, aerodynamic derivatives are all estimated with a central difference scheme, which will be presented later.Triangular membership functions.
(a)
(b)
(c)
(d)


## 2.3. Fuzzy Rule Inference

A fuzzy cell is formed by taking one membership function from each variable. The total number of cells is the number of possible combinations by taking one membership function from each input variable. For every cell, it has a fuzzy rule to guide the input and output relations. The rule of theith cell [15, 16] is stated as
(5)ifx1isA1i(x1),andx2isA2i(x2),and…,andxkisAki(xk),thenthecelloutputisPi=p0i+p1ix1+⋯+prixr+⋯+pkixk,
where i=1,2,…,n is the index of the cells, n is the total number of cells of the model; Pi(x1,x2,…,xr,…,xk) is the internal function with parameters p0i,p1i,…,pri,…,pki to be determined, and Aki(xk) denotes the membership function for xk. Each function covers a certain range of input variables.


## 2.4. Defuzzification

In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function, (5). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm (i.e., fuzzy rule inference). Because of this weighting among many factors over large ranges of possibilities, the word “fuzzy” is derived to describe the method. However, its prediction after defuzzification is in a deterministic value (a ‘‘crisp’’ one). The output estimated by the fuzzy-logic algorithm corresponding to the jth input (x1,j,x2,j,…,xr,j,…,xk,j) is as follows:
(6)y^j=∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]pi∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)].
In (6) product [Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)] is the weighted value of the ith 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 “product” stands product operator of its elements in this paper.


## 2.5. Parameter Identification

Given a set of membership functions for each input variable, the unknown coefficients of the internal functions are determined by using the gradient-descent method. The gradient descent is also known as steepest descent. The accuracy of the established aerodynamic model through the fuzzy-logic algorithm is estimated by the sum of squared errors (SSEs) and the multiple correlation coefficients(R2):
(7)SSE=∑j=1m(y^j-yj)2,(8)R2=1-{∑j=1m(y^j-yj)2}{∑j=1m(y--yj)2}.
In (7) and (8), where y^j, the output of the fuzzy-logic model at point j, is estimated by (6); yj is the data point used for the model training at point j; y- is the mean of the sample data, and m is the total number of data points. The model training is to determine the unknown coefficients of the internal functions, pri, by maximizing the value of R2. These coefficients are determined by the following iterative formula to minimize the sum of squared error (see (7)):
(9)pr,t+1i=pr,ti-αr∂(SSE)∂pri,
where αr is the convergence factor or the step size in the gradient method; subscript index t denotes the iteration sequence.After simplification, (9) becomes the following.Forr=0,(10a)p0,t+1i=p0,ti-2α0(y^j-yj)=×product[A1i(x1),…,Ari(xr,j),…,Aki(xk,j)]∑s=1nproduct[A1s(x1,j),…,Ars(xr,j),…,Aks(xk,j)],
and for r=1,2,…,k,
(10b)pr,t+1i=pr,ti-2αr(y^j-yj)=×product[A1i(x1),…,Ari(xr,j),…,Aki(xk,j)]xr,j∑s=1nproduct[A1s(x1,j),…,Ars(xr,j),…,Aks(xk,j)].The present FLM technique was utilized to conduct NASA research projects [17, 18]. The unsteady aerodynamic models for F-16XL were based on dynamic wind-tunnel data in [17], and those were from flight test data in [18]. The longitudinal and lateral-directional derivatives would be useful in aircraft stability analysis and control system design. Although the data source or application was not the same, the iteration during the search sequence stops of modeling should be the same with a principle for judgment.The iteration during the search sequence stops when one of the following three criteria is satisfied:(11)(1)Cost=SSEt<ε1,(2)RER=SSEt-SSEt-1SSEt<ε2,(3)t=tmax.
In the previos criteria, Cost=SSEI is the sum of squared errors (SSE) in current iteration to be denoted by “Cost” and RER=(cost_current-cost_previous)/cost_current to be denoted by “RER” (i.e., the relative error) for simplicity in descriptions; ε1 and ε2 are the required precision criteria; and tmax is a specified maximum iteration numbers.Given membership functions and the training data, this parameter identification procedure can be applied to establish a fuzzy-logic model. Although the procedure is well understood, to obtain a good model with the appropriate membership functions is quite challenging and requires numerical experimentation.


## 2.6. Model Structure Identification

In the fuzzy-logic model, the model structure is indicated by the number of membership functions for each variable. For a fuzzy-logic model with multiple variables, the structure is the combination of the numbers and forms of the membership functions assigned to all input variables. Since the sequence defines the one-to-one relationship between the numbers and the forms for each variable, the structure can be uniquely described by numbers of the membership functions.The model structure is determined by maximizing the correlation coefficient, (8). A search forward algorithm has been employed for the identification. At each search stage, there may be many fuzzy-logic models with different structure combinations. The search stage numbers are denoted by Ns. Out of all the possible intermediate fuzzy-logic models at each search stage, for an efficient search, only some structures are developed and evaluated. Two selection criteria, to be given below, are used to choose these structures. With the incremental sequence and the selection criteria, the search forward algorithm is summarized as follows.(1)
Specify the input variablesxr, r=1,2,…,k and the output variable y.(2)
Assume an initial structure also called parent structure as(N10,N20,…,Nr0,…,Nk0).(3)
Begin at the search stage numberNs=1, form all possible structures starting from the parent structure by adding one more membership function a time only to one input variable. Those all possible structures are called child structures as (N10+1,N20,…,Nr0,…,Nk0), (N10,N20+1,…,Nr0,…,Nk0),…,(N10,N20,…,Nr0,…,Nk0+1). Perform the identification of internal coefficients in (1) for each child structure and then calculate the R2 by using (8).(4)
Select the top 5 child structures among all calculated values ofR2 as new parent structures for next search step Ns=Ns+1.(5)
Go back to step(2) starting from the new parent structures and repeat the same procedures in steps (2) and (3) until the best structure is identified.(6)
Pick out the maximum value ofR2 among the child structures in each searching stage as Rmax2. The structure with the largest Rmax2 corresponding to all picking values is the optimal structure within a sensible Ns.The previous process is illustrated in Figure2. The structure identification is tied up with parameter identification. In the structure identification, parameter identification to determine the p-parameters according to (9) is also needed. The value of Nr is gradually changed in each iteration until the composite of Nr (i.e., optimal model structure) being fund, while the searching stage reaches Rmax2, for the model structure identification; and the composite of Nr is fixed, the number of iteration is unfixed until both the values of R2 and RER reach the requirements in the final the parameter identification.Figure 2
Identification process for the best structure.


## 3. Application to Aircraft Aerodynamic Modeling

### 3.1. Flight Data

In the present paper, a commercial transport aircraft encountered sever atmospheric turbulence twice during descending phase. The first and the second atmospheric turbulence encounters were at altitude about 10342 m and 9046 m, respectively. The dataset used for the modeling are extracted from the FDR during turbulence encounter lasting for 260 seconds.The main aircraft geometric and inertial characteristics for this aircraft are taken to be(12)W(take-off)=1,767,740N(397271lb),S=361.3m2(3890ft2),c-=6.005m(19.7ft),b=60.289m(197.8ft),Ixx=12,815,012kg-m2(9,452,600slugs-ft2),Iyy=12,511,197kg-m2(9,228,500slugs-ft2),Izz=53,650,986kg-m2(39,574,000slugs-ft2),Ixz=0.0kg-m2.The necessary data in the FDR to determine the aerodynamics for this aircraft is time(t), CAS, pressure altitude (h), Euler angles (ϕ, θ, and ψ), the longitudinal, lateral, and vertical accelerations (ax,ay,az), angle of attack (α), aileron deflection (δa), elevator δe, rudder δr, stabilizer δs, engine EPR, outside air temperature, wind speed (Vw), wind direction, and fuel flow rate. Since only the normal acceleration is recorded in 8-Hz resolution (i.e., 8 points per second), all other parameters are interpolated with a monotone cubic spline to the same sampling rate.


### 3.2. Compatibility Analysis

Typically, the longitudinal, lateral, and vertical accelerations(ax,ay,az) along the (x,y,z)-body axes of aircraft, angle of attack α, and the Euler angles (ϕ, θ, and ψ), as well as all control deflections are available and recorded in the FDR of all transport aircraft. Since the recorded flight data may contain errors (or called biases), compatibility analysis is performed to remove them by satisfying the following kinematic equations [19]:
(13)ϕ˙=p+qsinϕtanθ+rcosϕtanθ,θ˙=qcosϕ-rsinϕ,ψ˙=(qsinϕ+rcosϕ)secθ,V˙=(ax-gsinθ)cosαcosβ+(ay+gsinϕcosθ)sinβ+(az+gcosϕcosθ)sinαcosβ,α˙=[(az+gcosθcosϕ)cosα-(ax-gsinθ)sinα](Vcosβ)+q-tanβ(pcosα+rsinα),β˙=cosβ(ay+gcosθsinϕ)V+psinα-rcosα-sinβ[(az+gcosθcosϕ)sinα-(ax-gsinθ)cosα]V,
where g is the gravitational acceleration and V is the flight speed, β is sideslip angle, p is roll rate, q is pitch rate, and r is yaw rate in (13). Bias analysis is based on the least-square concept. The biases are: bax, bay, baz, bp, bq, br, bV, bα, bβ, bθ, bϕ, bψ, respectively, for ax, ay, az, and so forth, representing random noise and systematic errors, and are estimated by minimizing the sum of squares of differences between the two sides of the kinematic equations through the following algorithm:
(14)z⃑˙=f⃑(x⃑)=f→(x⃑m-Δx⃑),
where
(15)z⃑=(V,α-,β-,θ,ϕ,ψ)T,x⃑m=(a-x,a-y,a-z,p-,q-,r-,V,α-,β-,θ,ϕ,ψ)T,Δx⃑=(bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,bψ)T,
where super fix ‘‘-’’ stands for the mean value and the subscript “m” indicates the measured or recorded values. The cost function is defined as
(16)J=12(z⃑˙-f⃑)TQ(z⃑˙-f⃑),
where Q is a weighting diagonal matrix with elements being 1.0 except the one for the slowly varying flight speed being 10.0 and z⃑˙ is calculated with a central difference scheme with z⃑m, which is the measured value of z⃑. The steepest descent optimization method is adopted to minimize the cost function. As a result of the analysis, time derivatives, such as α˙, ϕ˙, and θ˙, are all estimated by using the central difference method; variables not present in the FDR, such as β, p, q, r, can be determined by cubic spline interpolation.Note that 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. Therefore the latter are obtained from the following flight dynamic equations [30]:
(17)max=Cxq-S+Tx,may=Cyq-S,maz=Czq-S,Clq-Sb=Ixxp˙-Ixzr˙+qr(Izz-Iyy)-Ixzpq,Cmq-Sc-=Iyyq˙+rp(Ixx-Izz)+Ixz(p2-r2)-Tm,Cnq-Sb=-Ixzp˙+Izzr˙+pq(Iyy-Ixx)+Ixzqr,
where m is the aircraft mass; q- is the dynamic pressure; S is the wing reference area; Cx, Cz, and Cm are the longitudinal aerodynamic force and moment coefficients; Cy, Cl, and Cn are the lateral-directional aerodynamic force and moment coefficients; and Ixx, Iyy, and Izz are the moments of inertia about x-, y-, and z-axes, respectively. The products of inertia, Ixy, Ixz, and Iyz, are assumed zero in the present case; but they are included in the equations because nonzero values may be available in other applications. The terms, Tx and Tm, represent the thrust contributions to the force in the direction of x-axes, and to the pitching moment, respectively.


### 3.3. Equivalent Harmonic Motion

The reduced frequency is a parameter to indicate the degree of unsteadiness in unsteady aerodynamics and is estimated in this paper by fitting the local trajectory with a harmonic motion. In the static case, the reduced frequency approaches 0. Large values of the reduced frequency imply the importance of unsteady aerodynamic effect. For longitudinal aerodynamics, the equivalent harmonic motion is the one based on the angle of attack variation (i.e., Theodorsen’s theory). For lateral-directional aerodynamics, it is based on the time variation of roll angle [20].For the longitudinal motion, the time history of the angle of attack(α) and time rate of angle of attack (dα/dt, or α˙) are fitted with one of a harmonic motion at any instant as follows:
(18)α(t)=α-+acos(ω1t+ϕ-),α˙(t)=-aω1sin(ω1t+ϕ-),
where those terms on the left hand side of (18) are given and the unknowns are the local mean angle of attack (α-), the local amplitude of the harmonic motion (a), the phase lag (ϕ-), and the estimated harmonic frequency (ω1). These unknowns are calculated through an optimization method by minimizing the following cost function (least squares):
(19)J=∑i=1n[αi-(α-+acos(ω1ti+ϕ-))]2+[α˙i-(α-+aω1sin(ω1ti+ϕ-))]2.
In (19), where αi is the measured value at point i and n is the number of the data points used in the optimization. For the case in the present study, n=20 is found to be the best value. The 20 points preceding and including the current time are employed in (19). The least-square method is found to converge well and gives reasonably accurate results. The lateral-directional equivalent reduced frequency is computed in the same manner.The local equivalent reduced frequency in the longitudinal motion is defined as(20)k1=ω1c-V,
where c- is the mean chord length of wing airfoil section; ω1 is the estimated harmonic frequency with time history of angle of attack.The lateral-directional equivalent reduced frequency is defined as(21)k2=ω2b2V,
where b is the wing span; ω2 is the estimated harmonic frequency with time history of roll angle.


### 3.4. Fuzzy-Logic Aerodynamic Models

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 atmospheric disturbances [19, 20]. For longitudinal aerodynamics, the models are assumed to be of the form
(22)Cx,Cz,Cm=f(α,α˙,q,k1,β,δe,M,p,δs,q-),
where the left hand side represents the coefficients of axial force (Cx), normal force (Cz), and pitching moment (Cm), respectively. All variables on the right hand side of (22) have been defined in the previous section. It should be noted that the stabilizer angle (δs) is included here, because it varies, though slowly, in flight to provide pitch trim (i.e., reducing the total static pitching moment to 0.0). The roll rate is included here because it is known that an aircraft under high aerodynamic loads at transonic speeds may have its longitudinal stability derivatives affected when additional disturbance due to roll rate is imposed.For the lateral-directional aerodynamics,(23)Cy,Cl,Cn=f(α,β,ϕ,p,r,k2,δa,δr,M,α˙,β˙),
where the left hand side represents the coefficients of side force (Cy), rolling moment (Cl), and yawing moment (Cn), respectively.


## 3.1. Flight Data

In the present paper, a commercial transport aircraft encountered sever atmospheric turbulence twice during descending phase. The first and the second atmospheric turbulence encounters were at altitude about 10342 m and 9046 m, respectively. The dataset used for the modeling are extracted from the FDR during turbulence encounter lasting for 260 seconds.The main aircraft geometric and inertial characteristics for this aircraft are taken to be(12)W(take-off)=1,767,740N(397271lb),S=361.3m2(3890ft2),c-=6.005m(19.7ft),b=60.289m(197.8ft),Ixx=12,815,012kg-m2(9,452,600slugs-ft2),Iyy=12,511,197kg-m2(9,228,500slugs-ft2),Izz=53,650,986kg-m2(39,574,000slugs-ft2),Ixz=0.0kg-m2.The necessary data in the FDR to determine the aerodynamics for this aircraft is time(t), CAS, pressure altitude (h), Euler angles (ϕ, θ, and ψ), the longitudinal, lateral, and vertical accelerations (ax,ay,az), angle of attack (α), aileron deflection (δa), elevator δe, rudder δr, stabilizer δs, engine EPR, outside air temperature, wind speed (Vw), wind direction, and fuel flow rate. Since only the normal acceleration is recorded in 8-Hz resolution (i.e., 8 points per second), all other parameters are interpolated with a monotone cubic spline to the same sampling rate.


## 3.2. Compatibility Analysis

Typically, the longitudinal, lateral, and vertical accelerations(ax,ay,az) along the (x,y,z)-body axes of aircraft, angle of attack α, and the Euler angles (ϕ, θ, and ψ), as well as all control deflections are available and recorded in the FDR of all transport aircraft. Since the recorded flight data may contain errors (or called biases), compatibility analysis is performed to remove them by satisfying the following kinematic equations [19]:
(13)ϕ˙=p+qsinϕtanθ+rcosϕtanθ,θ˙=qcosϕ-rsinϕ,ψ˙=(qsinϕ+rcosϕ)secθ,V˙=(ax-gsinθ)cosαcosβ+(ay+gsinϕcosθ)sinβ+(az+gcosϕcosθ)sinαcosβ,α˙=[(az+gcosθcosϕ)cosα-(ax-gsinθ)sinα](Vcosβ)+q-tanβ(pcosα+rsinα),β˙=cosβ(ay+gcosθsinϕ)V+psinα-rcosα-sinβ[(az+gcosθcosϕ)sinα-(ax-gsinθ)cosα]V,
where g is the gravitational acceleration and V is the flight speed, β is sideslip angle, p is roll rate, q is pitch rate, and r is yaw rate in (13). Bias analysis is based on the least-square concept. The biases are: bax, bay, baz, bp, bq, br, bV, bα, bβ, bθ, bϕ, bψ, respectively, for ax, ay, az, and so forth, representing random noise and systematic errors, and are estimated by minimizing the sum of squares of differences between the two sides of the kinematic equations through the following algorithm:
(14)z⃑˙=f⃑(x⃑)=f→(x⃑m-Δx⃑),
where
(15)z⃑=(V,α-,β-,θ,ϕ,ψ)T,x⃑m=(a-x,a-y,a-z,p-,q-,r-,V,α-,β-,θ,ϕ,ψ)T,Δx⃑=(bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,bψ)T,
where super fix ‘‘-’’ stands for the mean value and the subscript “m” indicates the measured or recorded values. The cost function is defined as
(16)J=12(z⃑˙-f⃑)TQ(z⃑˙-f⃑),
where Q is a weighting diagonal matrix with elements being 1.0 except the one for the slowly varying flight speed being 10.0 and z⃑˙ is calculated with a central difference scheme with z⃑m, which is the measured value of z⃑. The steepest descent optimization method is adopted to minimize the cost function. As a result of the analysis, time derivatives, such as α˙, ϕ˙, and θ˙, are all estimated by using the central difference method; variables not present in the FDR, such as β, p, q, r, can be determined by cubic spline interpolation.Note that 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. Therefore the latter are obtained from the following flight dynamic equations [30]:
(17)max=Cxq-S+Tx,may=Cyq-S,maz=Czq-S,Clq-Sb=Ixxp˙-Ixzr˙+qr(Izz-Iyy)-Ixzpq,Cmq-Sc-=Iyyq˙+rp(Ixx-Izz)+Ixz(p2-r2)-Tm,Cnq-Sb=-Ixzp˙+Izzr˙+pq(Iyy-Ixx)+Ixzqr,
where m is the aircraft mass; q- is the dynamic pressure; S is the wing reference area; Cx, Cz, and Cm are the longitudinal aerodynamic force and moment coefficients; Cy, Cl, and Cn are the lateral-directional aerodynamic force and moment coefficients; and Ixx, Iyy, and Izz are the moments of inertia about x-, y-, and z-axes, respectively. The products of inertia, Ixy, Ixz, and Iyz, are assumed zero in the present case; but they are included in the equations because nonzero values may be available in other applications. The terms, Tx and Tm, represent the thrust contributions to the force in the direction of x-axes, and to the pitching moment, respectively.


## 3.3. Equivalent Harmonic Motion

The reduced frequency is a parameter to indicate the degree of unsteadiness in unsteady aerodynamics and is estimated in this paper by fitting the local trajectory with a harmonic motion. In the static case, the reduced frequency approaches 0. Large values of the reduced frequency imply the importance of unsteady aerodynamic effect. For longitudinal aerodynamics, the equivalent harmonic motion is the one based on the angle of attack variation (i.e., Theodorsen’s theory). For lateral-directional aerodynamics, it is based on the time variation of roll angle [20].For the longitudinal motion, the time history of the angle of attack(α) and time rate of angle of attack (dα/dt, or α˙) are fitted with one of a harmonic motion at any instant as follows:
(18)α(t)=α-+acos(ω1t+ϕ-),α˙(t)=-aω1sin(ω1t+ϕ-),
where those terms on the left hand side of (18) are given and the unknowns are the local mean angle of attack (α-), the local amplitude of the harmonic motion (a), the phase lag (ϕ-), and the estimated harmonic frequency (ω1). These unknowns are calculated through an optimization method by minimizing the following cost function (least squares):
(19)J=∑i=1n[αi-(α-+acos(ω1ti+ϕ-))]2+[α˙i-(α-+aω1sin(ω1ti+ϕ-))]2.
In (19), where αi is the measured value at point i and n is the number of the data points used in the optimization. For the case in the present study, n=20 is found to be the best value. The 20 points preceding and including the current time are employed in (19). The least-square method is found to converge well and gives reasonably accurate results. The lateral-directional equivalent reduced frequency is computed in the same manner.The local equivalent reduced frequency in the longitudinal motion is defined as(20)k1=ω1c-V,
where c- is the mean chord length of wing airfoil section; ω1 is the estimated harmonic frequency with time history of angle of attack.The lateral-directional equivalent reduced frequency is defined as(21)k2=ω2b2V,
where b is the wing span; ω2 is the estimated harmonic frequency with time history of roll angle.


## 3.4. Fuzzy-Logic Aerodynamic Models

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 atmospheric disturbances [19, 20]. For longitudinal aerodynamics, the models are assumed to be of the form
(22)Cx,Cz,Cm=f(α,α˙,q,k1,β,δe,M,p,δs,q-),
where the left hand side represents the coefficients of axial force (Cx), normal force (Cz), and pitching moment (Cm), respectively. All variables on the right hand side of (22) have been defined in the previous section. It should be noted that the stabilizer angle (δs) is included here, because it varies, though slowly, in flight to provide pitch trim (i.e., reducing the total static pitching moment to 0.0). The roll rate is included here because it is known that an aircraft under high aerodynamic loads at transonic speeds may have its longitudinal stability derivatives affected when additional disturbance due to roll rate is imposed.For the lateral-directional aerodynamics,(23)Cy,Cl,Cn=f(α,β,ϕ,p,r,k2,δa,δr,M,α˙,β˙),
where the left hand side represents the coefficients of side force (Cy), rolling moment (Cl), and yawing moment (Cn), respectively.


## 4. Numerical Results and Discussions

The commercial transport aircraft encountered sever atmospheric turbulences in revenue flights in descending phase. To examine the static and dynamic stability characteristics, it is imperative to understand the flight environment in detail.


### 4.1. Aerodynamic and Flight Environments

The corresponding flight data for this aircraft in severe atmospheric turbulence during descending flight is presented in Figure3. The dataset of time span from t=7480~7739 sec used for the modeling is extracted from the FDR. This aircraft encounters severe atmospheric turbulence twice during this time span. In Figure 3(a), the variations of normal acceleration (az) show the highest az being 2.05 g around t=7483 sec and the lowest being −1.05 g around t=7484 sec in the first turbulence encounter, the highest az being 1.91 g around t=7682 sec and the lowest being −0.16 g around t=7684 sec in the second one. Figure 3(b) shows that variation of α is approximately in phase with az during those two turbulence encounters; α is highest about 4 deg. in the first turbulence encounter and 5.5 deg. in the second one. The altitude (h) with drop-off heights in the first turbulence encounter with the time span from t=7482~7494 sec is presented in Figure 3(c) and in the second turbulence encounter with the time span from t=7680~7690 sec is presented in Figure 3(d). The aircraft rapidly plunges downward during the turbulence encounter. The largest drop-off height reaches 57.3 m in the time span between t=7484~7486 sec, as shown in Figure 3(c). The Mach number (M) drops from 0.83 to 0.75 in the first turbulence encounter and from 0.80 to 0.70 in the second one, as shown in Figure 3(e).The time history of flight variables in severe atmospheric turbulence during the descending phase in transonic flight.
(a)
(b)
(c)
(d)
(e)The dynamic aerodynamic effects can be expected to be very significant under the circumstances of instantaneous changes ofα, h, and M in transonic flight. Since α reaches the value about 5.0 deg. at transonic flight, compressibility effect is important. It should be noted that the turbulent vertical wind field was not measured or estimated in the FDR; but it is included in the total α.


### 4.2. Pilot’s Responses in Operations

Figure4 presents the variations of dynamic characteristics and control variables. The variation ranges of α are 4 deg. to −1.8 deg. in the first turbulence encounter and 5.5 deg. to −2 deg. in the second one; the time history of β is about −2 deg. with the magnitude of small fluctuation, as indicated in Figure 4(a). The time history of pitch angle (θ) and roll angle (ϕ) is shown in Figure 4(b). The pitch angle (θ) does not vary as much as the angle of attack (α), reaching only 2.5 deg. The variations of ϕ are large, especially during the second turbulence encounter. The magnitudes of pitch rate (q) and roll rate (p) are shown in Figure 4(c); the variation of p is larger than that of q. The yaw rate is not shown because it is small throughout.The variations of dynamic characteristics and control variables in severe atmospheric turbulence during the descending phase in transonic flight.
(a)
(b)
(c)
(d)The flight data recorder indicated that the autopilot was deactivated fromt=7489 sec to t=7493 sec, for 4 seconds. It is of interest to examine the pilot’s response in operations after the autopilot being tripped off. The time history of δe and δa is shown in Figure 4(d). The magnitude of δe becomes positive to keep pitch down in descending. The variation of δa in the second turbulence encounter is larger than that of the first one. While the roll angle (ϕ) is changing fast in the negative direction right before t=7643 sec, δa is increasingly more positive, reaching +13 deg at t=7645 sec, as shown in Figure 4(d). On the other hand, opposite changes occur around t=7646 sec. The δr is not shown because it is small throughout.Both the pitch and roll rates are significant in severe atmospheric turbulence, especially the latter. The yaw rate is not shown in Figure4, because it is small throughout. In normal flight conditions, pitch and roll rates are typically small, in particular at transonic speeds, to avoid nonlinear aerodynamic effects. The autopilot was tripped off in severe atmospheric turbulence. The pitch angle and roll angle are the only cues to the pilots in making the control decision.


### 4.3. Analysis of Model Predictions

In the present study, the accuracy of the established unsteady aerodynamic models with six aerodynamic coefficients by using FLM technique is estimated by the sum of squared errors (SSE) and the square of multiple correlation coefficients(R2). All the aerodynamic derivatives in the study of stability and controllability are calculated with these aerodynamic models of aerodynamic coefficients.Figure5 presents the aerodynamic coefficients of normal force Cz, pitching moment Cm, rolling moment Cl, and yawing moment Cn predicted by the unsteady aerodynamic models. The predicted data by the final models has a good match with the flight data.The main aerodynamic coefficients predicted by the nonlinear and unsteady aerodynamic models.
(a)
(b)
(c)
(d)The final main aerodynamic models of aerodynamic coefficients consist of many fuzzy rules for each coefficient as described in Tables1 and 2. In Tables 1 and 2, the numbers below each input variable represents the number of membership function. The total number of fuzzy cells (n) in each model is the product of each number which presented in column 3. The last column shows the final multiple correlation coefficients (R2). The accuracy of the established aerodynamic model through the fuzzy-logic algorithm can be judged by the multiple correlation coefficients (R2).Table 1
The final main aerodynamic models of longitudinal aerodynamics.
Coef.
α
α
˙
q
k
1
β
δ
e
M
p
q
-
n
R
2
C
z
2
3
3
2
3
3
2
3
2
3888
0.9618
C
m
2
2
3
4
2
3
2
2
3
3456
0.9873Table 2
The final main aerodynamic models of lateral-directional aerodynamics.
Coef.
α
β
ϕ
p
r
k
2
δ
a
δ
r
M
α
˙
β
˙
n
R
2
C
l
2
3
3
2
2
3
2
2
3
3
2
15552
0.9617
C
n
3
2
4
2
2
2
3
3
2
2
2
13824
0.9435


### 4.4. Analysis of Stability Characteristics

Since the sudden plunging motion with the abrupt change in altitude affects the flight safety the most in severe atmospheric turbulence, so the time period of7482~7494 sec is emphasized in examining the characteristics of stability. In order to examine the variations in characteristics, the units of all aerodynamic derivatives are converted to rad−1. It should be noted that these derivatives are evaluated at the instantaneous conditions, instead of about the trim conditions as have been traditionally done.The time history of main longitudinal and lateral-directional static stability derivatives along the flight path is presented in Figure6. From the point of view in static stability, initially, the configuration has longitudinal stability (Czα>0 and Cmα<0) in Figure 6(a), stable longitudinal damping (Cmq<0) in Figure 6(b), lateral stability (Clβ<0) and directional stability (Cnβ>0) in Figure 6(c), roll damping (Clp<0), and directional damping (Cnr<0) in Figure 6(d). During the plunging motion, Cmα>0 and Cmq>0 in the period of t=7483.0~7484.5 sec, and the static stability becomes unstable with insufficient roll damping (Clp≈0) and directional damping (Cnr>0). The aerodynamic instability is most likely caused by the motion that produces a time-dependent pressure distribution on the aircraft surface involving compressibility effects.The time history of main longitudinal and lateral-directional static stability derivatives along the flight path.
(a)
(b)
(c)
(d)Figure7 presents the time history of main longitudinal and lateral-directional oscillatory derivatives along the flight path to associate with α˙ and β˙-derivatives. Note that in Figure 7(a), the oscillatory derivatives are defined as
(24)(Cmq)osc=Cmq+Cmα˙,(Czq)osc=Czq+Czα˙.The time history of main longitudinal and lateral-directional oscillatory derivatives along the flight path.
(a)
(b)
(c)
(d)In Figure7(c), the oscillatory derivatives are defined as
(25)(Clp)osc=Clp+Clβ˙sinα,(Cnr)osc=Cnr-Cnβ˙cosα.The values of oscillatory derivatives are equivalent to the combinations of steady damping and dynamic derivatives in the previous equations (24)~(25). The use of oscillatory derivative instead of steady damping ones is more consistent with the actual case of the aerodynamic damping in the analysis of stability characteristics. To be stable, (Czq)osc>0, (Cmq)osc<0, (Clp)osc<0, and (Cnr)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 steady damping derivatives in Figures 7(a) and 7(c) due to the effects of the dynamic derivatives (i.e., α˙- and β˙-derivatives).Figures7(b) and 7(d) show dynamic derivatives of stability. To be stable, Czα˙>0, Cmα˙<0, Clβ˙<0, and Cnβ˙>0. The magnitudes of Czα˙ and Cmα˙ have significant variations and Czα˙<0 in the period of t=7482.5~7491 sec in Figure 7(b). It should be noted that Czα˙ represents the virtual mass effect and is particularly large in transonic flow to affect the plunging motion [31]. The (Czq)osc and (Cmq)osc are inadequate in oscillatory damping in the periods of t=7484~7486 sec and t=7490~7494 sec, as shown in Figure 7(a). The effect of α˙-derivative on (Cmq)osc is to improve the stability in pitch in the periods of t=7482.5~7484 sec and 7486~7488.5 sec. The magnitudes of Clβ˙ in Figure 7(d) and (Clp)osc in Figure 7(c) are small; Clp is not shown because it is small throughout. The values of Cnβ˙ are from positive at t=7483 sec to negative at t=7484 sec; (Cnr)osc are from positive at t=7483.5 sec to negative at t=7485 sec. It implies that the effects of β˙-derivative are to cause the directional stability more unstable.In essence, the effects ofα˙-derivative on (Czq)osc and β˙-derivative on (Clp)osc are small. However, the effect of α˙-derivative on (Cmq)osc is to improve the stability in pitch; while the effects of β˙-derivative are to cause more directional instability. These results indicate that the turbulent crosswind has some adverse effects on directional stability and damping. Although the dynamic derivatives tend to be small for the present configuration, this is much helpful to understand the unknown factors of instability characteristics.


## 4.1. Aerodynamic and Flight Environments

The corresponding flight data for this aircraft in severe atmospheric turbulence during descending flight is presented in Figure3. The dataset of time span from t=7480~7739 sec used for the modeling is extracted from the FDR. This aircraft encounters severe atmospheric turbulence twice during this time span. In Figure 3(a), the variations of normal acceleration (az) show the highest az being 2.05 g around t=7483 sec and the lowest being −1.05 g around t=7484 sec in the first turbulence encounter, the highest az being 1.91 g around t=7682 sec and the lowest being −0.16 g around t=7684 sec in the second one. Figure 3(b) shows that variation of α is approximately in phase with az during those two turbulence encounters; α is highest about 4 deg. in the first turbulence encounter and 5.5 deg. in the second one. The altitude (h) with drop-off heights in the first turbulence encounter with the time span from t=7482~7494 sec is presented in Figure 3(c) and in the second turbulence encounter with the time span from t=7680~7690 sec is presented in Figure 3(d). The aircraft rapidly plunges downward during the turbulence encounter. The largest drop-off height reaches 57.3 m in the time span between t=7484~7486 sec, as shown in Figure 3(c). The Mach number (M) drops from 0.83 to 0.75 in the first turbulence encounter and from 0.80 to 0.70 in the second one, as shown in Figure 3(e).The time history of flight variables in severe atmospheric turbulence during the descending phase in transonic flight.
(a)
(b)
(c)
(d)
(e)The dynamic aerodynamic effects can be expected to be very significant under the circumstances of instantaneous changes ofα, h, and M in transonic flight. Since α reaches the value about 5.0 deg. at transonic flight, compressibility effect is important. It should be noted that the turbulent vertical wind field was not measured or estimated in the FDR; but it is included in the total α.


## 4.2. Pilot’s Responses in Operations

Figure4 presents the variations of dynamic characteristics and control variables. The variation ranges of α are 4 deg. to −1.8 deg. in the first turbulence encounter and 5.5 deg. to −2 deg. in the second one; the time history of β is about −2 deg. with the magnitude of small fluctuation, as indicated in Figure 4(a). The time history of pitch angle (θ) and roll angle (ϕ) is shown in Figure 4(b). The pitch angle (θ) does not vary as much as the angle of attack (α), reaching only 2.5 deg. The variations of ϕ are large, especially during the second turbulence encounter. The magnitudes of pitch rate (q) and roll rate (p) are shown in Figure 4(c); the variation of p is larger than that of q. The yaw rate is not shown because it is small throughout.The variations of dynamic characteristics and control variables in severe atmospheric turbulence during the descending phase in transonic flight.
(a)
(b)
(c)
(d)The flight data recorder indicated that the autopilot was deactivated fromt=7489 sec to t=7493 sec, for 4 seconds. It is of interest to examine the pilot’s response in operations after the autopilot being tripped off. The time history of δe and δa is shown in Figure 4(d). The magnitude of δe becomes positive to keep pitch down in descending. The variation of δa in the second turbulence encounter is larger than that of the first one. While the roll angle (ϕ) is changing fast in the negative direction right before t=7643 sec, δa is increasingly more positive, reaching +13 deg at t=7645 sec, as shown in Figure 4(d). On the other hand, opposite changes occur around t=7646 sec. The δr is not shown because it is small throughout.Both the pitch and roll rates are significant in severe atmospheric turbulence, especially the latter. The yaw rate is not shown in Figure4, because it is small throughout. In normal flight conditions, pitch and roll rates are typically small, in particular at transonic speeds, to avoid nonlinear aerodynamic effects. The autopilot was tripped off in severe atmospheric turbulence. The pitch angle and roll angle are the only cues to the pilots in making the control decision.


## 4.3. Analysis of Model Predictions

In the present study, the accuracy of the established unsteady aerodynamic models with six aerodynamic coefficients by using FLM technique is estimated by the sum of squared errors (SSE) and the square of multiple correlation coefficients(R2). All the aerodynamic derivatives in the study of stability and controllability are calculated with these aerodynamic models of aerodynamic coefficients.Figure5 presents the aerodynamic coefficients of normal force Cz, pitching moment Cm, rolling moment Cl, and yawing moment Cn predicted by the unsteady aerodynamic models. The predicted data by the final models has a good match with the flight data.The main aerodynamic coefficients predicted by the nonlinear and unsteady aerodynamic models.
(a)
(b)
(c)
(d)The final main aerodynamic models of aerodynamic coefficients consist of many fuzzy rules for each coefficient as described in Tables1 and 2. In Tables 1 and 2, the numbers below each input variable represents the number of membership function. The total number of fuzzy cells (n) in each model is the product of each number which presented in column 3. The last column shows the final multiple correlation coefficients (R2). The accuracy of the established aerodynamic model through the fuzzy-logic algorithm can be judged by the multiple correlation coefficients (R2).Table 1
The final main aerodynamic models of longitudinal aerodynamics.
Coef.
α
α
˙
q
k
1
β
δ
e
M
p
q
-
n
R
2
C
z
2
3
3
2
3
3
2
3
2
3888
0.9618
C
m
2
2
3
4
2
3
2
2
3
3456
0.9873Table 2
The final main aerodynamic models of lateral-directional aerodynamics.
Coef.
α
β
ϕ
p
r
k
2
δ
a
δ
r
M
α
˙
β
˙
n
R
2
C
l
2
3
3
2
2
3
2
2
3
3
2
15552
0.9617
C
n
3
2
4
2
2
2
3
3
2
2
2
13824
0.9435


## 4.4. Analysis of Stability Characteristics

Since the sudden plunging motion with the abrupt change in altitude affects the flight safety the most in severe atmospheric turbulence, so the time period of7482~7494 sec is emphasized in examining the characteristics of stability. In order to examine the variations in characteristics, the units of all aerodynamic derivatives are converted to rad−1. It should be noted that these derivatives are evaluated at the instantaneous conditions, instead of about the trim conditions as have been traditionally done.The time history of main longitudinal and lateral-directional static stability derivatives along the flight path is presented in Figure6. From the point of view in static stability, initially, the configuration has longitudinal stability (Czα>0 and Cmα<0) in Figure 6(a), stable longitudinal damping (Cmq<0) in Figure 6(b), lateral stability (Clβ<0) and directional stability (Cnβ>0) in Figure 6(c), roll damping (Clp<0), and directional damping (Cnr<0) in Figure 6(d). During the plunging motion, Cmα>0 and Cmq>0 in the period of t=7483.0~7484.5 sec, and the static stability becomes unstable with insufficient roll damping (Clp≈0) and directional damping (Cnr>0). The aerodynamic instability is most likely caused by the motion that produces a time-dependent pressure distribution on the aircraft surface involving compressibility effects.The time history of main longitudinal and lateral-directional static stability derivatives along the flight path.
(a)
(b)
(c)
(d)Figure7 presents the time history of main longitudinal and lateral-directional oscillatory derivatives along the flight path to associate with α˙ and β˙-derivatives. Note that in Figure 7(a), the oscillatory derivatives are defined as
(24)(Cmq)osc=Cmq+Cmα˙,(Czq)osc=Czq+Czα˙.The time history of main longitudinal and lateral-directional oscillatory derivatives along the flight path.
(a)
(b)
(c)
(d)In Figure7(c), the oscillatory derivatives are defined as
(25)(Clp)osc=Clp+Clβ˙sinα,(Cnr)osc=Cnr-Cnβ˙cosα.The values of oscillatory derivatives are equivalent to the combinations of steady damping and dynamic derivatives in the previous equations (24)~(25). The use of oscillatory derivative instead of steady damping ones is more consistent with the actual case of the aerodynamic damping in the analysis of stability characteristics. To be stable, (Czq)osc>0, (Cmq)osc<0, (Clp)osc<0, and (Cnr)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 steady damping derivatives in Figures 7(a) and 7(c) due to the effects of the dynamic derivatives (i.e., α˙- and β˙-derivatives).Figures7(b) and 7(d) show dynamic derivatives of stability. To be stable, Czα˙>0, Cmα˙<0, Clβ˙<0, and Cnβ˙>0. The magnitudes of Czα˙ and Cmα˙ have significant variations and Czα˙<0 in the period of t=7482.5~7491 sec in Figure 7(b). It should be noted that Czα˙ represents the virtual mass effect and is particularly large in transonic flow to affect the plunging motion [31]. The (Czq)osc and (Cmq)osc are inadequate in oscillatory damping in the periods of t=7484~7486 sec and t=7490~7494 sec, as shown in Figure 7(a). The effect of α˙-derivative on (Cmq)osc is to improve the stability in pitch in the periods of t=7482.5~7484 sec and 7486~7488.5 sec. The magnitudes of Clβ˙ in Figure 7(d) and (Clp)osc in Figure 7(c) are small; Clp is not shown because it is small throughout. The values of Cnβ˙ are from positive at t=7483 sec to negative at t=7484 sec; (Cnr)osc are from positive at t=7483.5 sec to negative at t=7485 sec. It implies that the effects of β˙-derivative are to cause the directional stability more unstable.In essence, the effects ofα˙-derivative on (Czq)osc and β˙-derivative on (Clp)osc are small. However, the effect of α˙-derivative on (Cmq)osc is to improve the stability in pitch; while the effects of β˙-derivative are to cause more directional instability. These results indicate that the turbulent crosswind has some adverse effects on directional stability and damping. Although the dynamic derivatives tend to be small for the present configuration, this is much helpful to understand the unknown factors of instability characteristics.


## 5. Concluding Remarks

The main objective in this paper was to present the development of nonlinear and unsteady aerodynamic models based on FLM method and the resulting models having the capability to create the required continuous derivatives for the investigative analysis. The static and dynamic stability characteristics in severe atmospheric turbulence with sudden plunging motion in descending flight were examined by using those continuous derivatives. The numerical results and discussions of application to aerodynamic modeling were concluded as follows.(1)
The uncertainty and correlation of the data points were estimated and improved by monitoring a multivariable correlation coefficient in the modeling process.(2)
The predicted results showed that the final models could not only produce relatively accurate aerodynamic coefficients but also generate smooth coefficient curves with the required independent variables to provide reasonable local aerodynamic derivatives.(3)
The fuzzy logic-based aerodynamic modeling method was shown to be capable of handling nonlinear and unsteady aerodynamic environments possible exhibited in abnormal flight condition.The analytical results could provide the mitigation concepts and promote the understanding of aerodynamic responses of the commercial transport aircraft in severe atmospheric turbulence.



---

*Source: 609769-2013-06-06.xml*