# Nonlinear and Dynamic Aerodynamic Models  for Commercial Transport Aircraft with Adverse Weather Effects

Authors: Ray C. Chang; Cun-En Ye; C. Edward Lan; Yannian Lee
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The flight dynamic equations in mathematics for aircraft response to the adverse weathers, such as wind shear, atmospheric turbulence, and in-flight icing, are nonlinear and unsteady. To effectively analyze the performance degradation and variations in stability of commercial aircraft that encountered these weather hazards, the nonlinear and dynamic (i.e., time dependent) aerodynamic models based on flight data would be needed. In the present paper, a numerical modeling method based on a fuzzy-logic algorithm will be presented to estimate the aerodynamic models for a twin-jet transport by using the flight data from the flight data recorder (FDR). The aerodynamic models having the capability to generate continuous stability derivatives, especially for sensitivity study of unknown factors in adverse weather conditions, will be demonstrated in this paper.


---

## Body

## 1. Introduction

Aircrafts in flight are frequently subject to atmospheric disturbances. The hazards due to these disturbances may be in the form of wind shear, turbulence (both clear-air and convective), thunderstorms, in-flight icing, and so forth. The severity of aircraft response to the disturbance is related to the dynamic aerodynamics that results from the instantaneous changes of aircraft flight attitudes. To provide the mitigation concepts and promote the understanding of aerodynamic responses of the commercial transport aircraft in adverse weather conditions, the nonlinear and dynamic (i.e., time-dependent) aerodynamic models based on flight data would be needed.Regarding the flight data of commercial transport aircraft, the flight data recorder (FDR) is ICAO-regulated devices. The FDR, popularly referred to as a “black box”, is a mandatory device used to record specific aircraft parameters for event investigations. Although the data resolution and accuracy in the FDR are not as well defined as those in flight test, the predicted aerodynamic derivatives correlated well with the nominal values of transport aircraft in normal flight [1]. Therefore, the FDR data will be used in the present analysis.In setting up the aerodynamic models to predict the required derivatives in aerodynamics, the traditional method of using the tabulated data derived from such methods as the maximum likelihood method (MMLE) [2], the least-square or the stepwise regression method [3], is very difficult to correlate complex functional relations among numerous parameters. In earlier development of the fuzzy-logic algorithm, Zadeh [4] 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 is not suitable when derivatives of the function are needed. In 1985, Takagi and Sugeno [5] used the internal functions, instead of the fuzzy sets, in developing the fuzzy-logic algorithm. In this approach, the main features are the internal functions, membership functions, and the output cells. The output curves in the prediction are smooth. The present development will be based on this second approach. This second approach, to be called the Fuzzy-Logic Modeling (FLM) technique [6, 7], is capable of correlating multiple parameters without assuming explicit functional relations.Among the different adverse weathers, clear-air turbulence is all the more important in flight safety since it is hard to detect and predict. Clear-air turbulence is the leading cause of serious personal injuries in nonfatal accidents of commercial aircraft. One main type of motion that causes flight injuries in clear-air turbulence is the sudden plunging motion with the abrupt change in altitude. This paper presents FLM technique to establish unsteady aerodynamic models with six aerodynamic coefficients based on the datasets from the flight data recorder (FDR) of a twin-jet transport. The aerodynamic derivatives extracted from these aerodynamic models are then used to evaluate the variations in stability. The sensitivity study of unknown factors during the sudden plunging motion in severe clear-air turbulence will be demonstrated in this paper.


## 2. Theoretical Development

The general idea of the 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 [6–8].


### 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 [8] as follows:
(2.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; i=1,2,…,n, and n is the total number of fuzzy cells.The recorded data in FDR, such as flight altitude(h), calibrated airspeed (CAS), angle of attack (α), accelerometer readings (ax, ay, and az), Euler angles (θ, ϕ, and ψ), and so forth, 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, with “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 grades. In this paper, the overlapped triangular membership function is used to represent the grades of internal functions.The membership functions partition the input space into many fuzzy subspaces, which are called the fuzzy cells. 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 let j be the index for the jth membership functions. Then the membership grades can be described as follows [9]:(1)
N=2, A(xr)=xr, j=1,
(2.2)A(xr)=1-xr,j=2,(2)
j=3 to m, where m=max(0,(N-2)/2):
(2.3)A(xr)=xrdu,0≤xr≤du,A(xr)=(1-xr)(1-du),du≤xr≤1,
where du=Δx1*(j-2), and Δx1=1.0/(N-m-1),(3)
j≥N-m,(2.4)A(xr)=(dd-xr)dd,0≤xr≤dd,A(xr)=(dd-xr)(dd-1),dd≤xr≤1,
where dd = Δx2*(j-N+m), and Δx2=1.0/(m+1).


### 2.3. Fuzzy Rule and Output

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 the ith cell is stated [8] as follows: ifx1isA1i(x1), if  x2  isA2i(x2),…, and ifxk  isAki(xk), then the cell output is(2.5)Pi(x1,x2,…,xr,…xk)=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 grade for xk. Each function covers a certain range of input variables.In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function (2.5). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm. Because of this weighting among many factors over large ranges of possibilities, the word “fuzzy” is derived to describe the method. However, its prediction is never “fuzzy”. The output estimated by the fuzzy-logic algorithm corresponding to the jth input (x1,j,x2,j,…,xr,j,…,xk,j) is as follows:(2.6)ŷj=∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]pi∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)].
In (2.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,m is the total number of the data record, and the “product” stands for product operator of its elements in this paper.


### 2.4. Parameter Identification

Given a set of membership functions for each input variable, the unknown coefficients of the internal functions are adjusted by using the Newton gradient-descent method. 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):(2.7)SSE=∑j=1m(ŷj-yj)2,(2.8)R2=1-{∑j=1m(ŷj-yj)2}{∑j=1m(y̅-yj)2}.In (2.7) and (2.8), where ŷj, the output of the fuzzy-logic model at point j, is estimated by (2.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 (2.7):
(2.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.In usual, a large value ofαr is chosen, such as 10, at the beginning to accelerate the process of convergence; the small value of less then 1 is taken to obtain more precise unknown coefficients of the internal functions, pri, in model training process. After simplification, (2.9) becomes as follows.Forr=0,(2.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,(2.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 iteration during the search sequence stops when one of the following three criteria [6, 9] is satisfied:(1)
Cost=SSEt<ε1,(2)
RER=(SSEt-SSEt-1)/SSEt<ε2,(3)
t=tmax.In the above criteria,Cost=SSEt is the sum of squared errors (SSEs) 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 number.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 that  are quite challenging, it  requires numerical experimentation.


### 2.5. Model Structure Identification

In the fuzzy-logic model, the model structure is indicated by the number of membership functions. 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 the numbers.The model structure is optimized by maximizing (2.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 can be shown graphically in Figure 1.Figure 1
Best model structure searching flow.The step-by-step searching flow is summarized as follows [6].(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=l 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),…, and (N10,N20,…,Nr0,…,Nk0+1). Do the identification of internal coefficients in (2.1) for each child structure and then calculate its R2 by using (2.8).(4)
Select the top 5 child structures among all calculated values ofR2 as new parent structures for next search step Ns=Ns+l.(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 of the largest Rmax2 corresponding to the all picking values is the optimal structure within a sensible Ns.In the structure identification, parameter identification to determine thep-parameters according to (2.9) is also needed. But the number of iteration to determine the p-parameters is limited to 2000, so that the best structure is decided on a relative basis. After this last step, (2.9) is applied iteratively until both values of R2 and RER reach the requirements in the final parameter identification


## 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 [8] as follows:
(2.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; i=1,2,…,n, and n is the total number of fuzzy cells.The recorded data in FDR, such as flight altitude(h), calibrated airspeed (CAS), angle of attack (α), accelerometer readings (ax, ay, and az), Euler angles (θ, ϕ, and ψ), and so forth, 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, with “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 grades. In this paper, the overlapped triangular membership function is used to represent the grades of internal functions.The membership functions partition the input space into many fuzzy subspaces, which are called the fuzzy cells. 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 let j be the index for the jth membership functions. Then the membership grades can be described as follows [9]:(1)
N=2, A(xr)=xr, j=1,
(2.2)A(xr)=1-xr,j=2,(2)
j=3 to m, where m=max(0,(N-2)/2):
(2.3)A(xr)=xrdu,0≤xr≤du,A(xr)=(1-xr)(1-du),du≤xr≤1,
where du=Δx1*(j-2), and Δx1=1.0/(N-m-1),(3)
j≥N-m,(2.4)A(xr)=(dd-xr)dd,0≤xr≤dd,A(xr)=(dd-xr)(dd-1),dd≤xr≤1,
where dd = Δx2*(j-N+m), and Δx2=1.0/(m+1).


## 2.3. Fuzzy Rule and Output

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 the ith cell is stated [8] as follows: ifx1isA1i(x1), if  x2  isA2i(x2),…, and ifxk  isAki(xk), then the cell output is(2.5)Pi(x1,x2,…,xr,…xk)=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 grade for xk. Each function covers a certain range of input variables.In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function (2.5). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm. Because of this weighting among many factors over large ranges of possibilities, the word “fuzzy” is derived to describe the method. However, its prediction is never “fuzzy”. The output estimated by the fuzzy-logic algorithm corresponding to the jth input (x1,j,x2,j,…,xr,j,…,xk,j) is as follows:(2.6)ŷj=∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]pi∑i=1nproduct[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)].
In (2.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,m is the total number of the data record, and the “product” stands for product operator of its elements in this paper.


## 2.4. Parameter Identification

Given a set of membership functions for each input variable, the unknown coefficients of the internal functions are adjusted by using the Newton gradient-descent method. 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):(2.7)SSE=∑j=1m(ŷj-yj)2,(2.8)R2=1-{∑j=1m(ŷj-yj)2}{∑j=1m(y̅-yj)2}.In (2.7) and (2.8), where ŷj, the output of the fuzzy-logic model at point j, is estimated by (2.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 (2.7):
(2.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.In usual, a large value ofαr is chosen, such as 10, at the beginning to accelerate the process of convergence; the small value of less then 1 is taken to obtain more precise unknown coefficients of the internal functions, pri, in model training process. After simplification, (2.9) becomes as follows.Forr=0,(2.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,(2.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 iteration during the search sequence stops when one of the following three criteria [6, 9] is satisfied:(1)
Cost=SSEt<ε1,(2)
RER=(SSEt-SSEt-1)/SSEt<ε2,(3)
t=tmax.In the above criteria,Cost=SSEt is the sum of squared errors (SSEs) 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 number.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 that  are quite challenging, it  requires numerical experimentation.


## 2.5. Model Structure Identification

In the fuzzy-logic model, the model structure is indicated by the number of membership functions. 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 the numbers.The model structure is optimized by maximizing (2.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 can be shown graphically in Figure 1.Figure 1
Best model structure searching flow.The step-by-step searching flow is summarized as follows [6].(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=l 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),…, and (N10,N20,…,Nr0,…,Nk0+1). Do the identification of internal coefficients in (2.1) for each child structure and then calculate its R2 by using (2.8).(4)
Select the top 5 child structures among all calculated values ofR2 as new parent structures for next search step Ns=Ns+l.(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 of the largest Rmax2 corresponding to the all picking values is the optimal structure within a sensible Ns.In the structure identification, parameter identification to determine thep-parameters according to (2.9) is also needed. But the number of iteration to determine the p-parameters is limited to 2000, so that the best structure is decided on a relative basis. After this last step, (2.9) is applied iteratively until both values of R2 and RER reach the requirements in the final parameter identification


## 3. Application to Aerodynamic Modeling

### 3.1. Flight Data

The twin-jet transport of the present study encountered clear-air turbulence in cruise flight at the altitude around 10,050 m. As a result, several passengers and cabin crews sustained injuries, because of which this event was classified as an accident. The present study was initiated to examine possible concepts of accident prevention in the future. The datasets used for the modeling are extracted from the FDR during turbulence encounter lasting for 92 seconds.The main aircraft geometric and inertial characteristics are taken to beW (take-off) = 1,431,800 N (321900 lb),S= 260 m2  (2798.7 ft2), c̅= 6.608 m (21.68 ft), and b= 44.827 m (147.08 ft),Ixx= 10,710,000 kg-m2 (7,899,900 slugs-ft2), and Iyy= 14,883,800 kg-m2 (10,978,000 slugs-ft2),Izz= 25,283,271 kg-m2 (18,648,470 slugs-ft2) and Ixz= 0.0 kg-m2.The required operational parameters in FDR dataset for generating aerodynamic data files are time (t), CAS, pressure altitude (h), roll attitude (ϕ), pitch attitude (θ), magnetic heading (ψ), longitudinal acceleration (ax), lateral acceleration (ay), vertical acceleration (az), angle of attack (α), aileron deflection (δa), elevator (δe), rudder (δr), stabilizer (δs), engine EPR, outside air temperature, wind speed, 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:(3.1)ϕ̇=p+qsinϕtanθ+rcosϕtanθ,(3.2)θ̇=qcosϕ-rsinϕ,(3.3)ψ̇=(qsinϕ+rcosϕ)secθ,(3.4)V̇=(ax-gsinθ)cosαcosβ+(ay+gsinϕcosθ)sinβ+(az+gcosϕcosθ)sinαcosβ,(3.5)α̇=[(az+gcosθcosϕ)cosα-(ax-gsinθ)sinα](Vcosβ)+q-tanβ(pcosα+rsinα),(3.6)β̇=cosβ(ay+gcosθsinϕ)V+psinα-rcosα-sinβ[(az+gcosθcosϕ)sinα-(ax-gsinθ)cosα]V,
where g is acceleration due to gravity, V is flight speed, β is sideslip angle, p is roll rate, q is pitch rate, and r is yaw rate in (3.1)–(3.6). Let the biases be denoted by bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,andbψ, respectively, for ax, ay, az, and so forth. These biases are estimated by minimizing the squared sum of the differences between the two sides of the above equations. These equations in vector form can be written as(3.7)z⇀̇=f⇀(x)=f⇀(xm-Δx),
where(3.8)z⇀=(V,α,β,θ,ϕ,ψ)T,(3.9)x⇀m=(ax,ay,az,p,q,r,V,α,β,θ,ϕ,ψ)T,(3.10)Δx⇀=(bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,bφ)T,
where the subscript “m” indicates the measured or recorded values. The cost function is defined as(3.11)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, variables not present in the FDR, such as β, p, q, and r, are also estimated.The force and moment coefficients are obtained from the following flight dynamic equations [7, 10] about the airplane body axes:(3.12)max=Cxq̅S+Tx,(3.13)may=Cyq̅S+Ty,(3.14)maz=Czq̅S+Tz,(3.15)Clq̅Sb=Ixxṗ-Ixz(ṙ+pq)-(Iyy-Izz)qr,(3.16)Cmq̅Sc̅=Iyyq̇-Ixz(r2-p2)-(Izz-Ixx)rp-Tm,(3.17)Cnq̅Sb=Izzṙ-Ixz(ṗ-qr)-(Ixx-Iyy)pq,
where m is denoted as aircraft mass; q̅ is denoted as dynamic pressure; S is denoted as wing reference area; Cx, Cz, and Cm are denoted as longitudinal aerodynamic force and moment coefficients; Cy, Cl, and Cn are denoted as lateral-directional aerodynamic force and moment coefficients; Ixx, Iyy, and Izz are denoted as moments of inertia about x-, y-, and z-axes, respectively; Ixy, Ixz, and Iyz are denoted as products of inertia; Tx, Ty, Tz, and Tm are denoted as thrust terms about x-, y-, z-axes, and in equation of pitching moment, respectively, in (3.7)–(3.12).The above equations are used to determine all aerodynamic coefficients based on accelerometer readings (ax, ay, and az), Euler angles (ϕ, θ, and ψ), angular rates (p, q, and r), and thrusts (Tx, Ty, Tz, and Tm). The angular rates are estimated through compatibility analysis. Since thrusts were not measured during flight for most flight vehicles, those values and the effects on the forces and pitching moments in equations of (3.7), (3.8), (3.9), and (3.11) should be predicted by a thrust model.


### 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 is 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 following the classical unsteady aerodynamic theory of Theodorsen [11]. For lateral-directional aerodynamics, it is based on the time variation of roll angle [8].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 [7, 12] as follows:(3.18)α(t)=α̅+αcos(ωt+ϕ̅),(3.19)α̇(t)=-aωsin(ωt+ϕ̅),
where those terms on the left hand side of (3.13) and (3.14) 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 angular frequency (ω). These unknowns are calculated through an optimization method by minimizing the following cost function (least squares):(3.20)J=∑i=1n[αi-(α̅+acos(ωti+ϕ̅))]2+[α̇i-(α̅+aωsin(ωti+ϕ̅))]2.In (3.15), α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 (3.15). 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(3.21)k1=ωc̅V,
where c̅ is the mean chord length of wing airfoil section.The lateral-directional equivalent reduced frequency is defined as(3.22)k2=ωbV,
where b is the wing span.


### 3.4. Fuzzy-Logic Thrust Model

As shown before, the thrust terms appear in the force equations and the pitching moment equations (3.7)–(3.9) and (3.11). Since the values of thrust for aircraft in flight cannot be directly measured in the current state of the art, they are not recorded in the FDR. The manufacturers of airplanes and engines and safety agencies agreed that using such parameters as the Mach number, airspeed, flight altitude, temperature, the rpm of the pressure compressors, and engine pressure ratios is adequate to estimate the engine thrust. 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 realistic one tied to the recorded engine performance parameters will be developed with the fuzzy-logic algorithm.For a commercial aircraft, most likely only the axial force and the pitching moment are affected by thrust. This assumption will be made in this paper. Theoretically, clear-air turbulence (i.e., random change inu, w (or α), and v (or β) affects the engine performance through its effects on static and dynamic distortions at the engine face. However, the effects of clear-air turbulence on the engine parameters and the thrust magnitude are not known and, therefore, are ignored.For this purpose, data from the flight manual for the fuel flow rates (ṁf) at various altitudes (h), weights (W), Mach numbers (M), calibrated airspeed (CAS), engine pressure ratios (EPRs) in cruise flight are utilized. Note that the drag polar for a given aircraft is generally not known to most researchers. To estimate it and, hence, the thrust magnitude in cruise, the assumption of a design lift-to-drag ratio (L/D) of 17.5 is made. This value of lift-to-drag in cruise is assumed based on the past design experience for twin-jet transports. In the flight manual, various weights, altitudes, Mach numbers, CAS, EPR, and fuel flow rates in cruise are tabulated. The lift coefficient can be calculated at each flight condition immediately. As a result, the drag coefficient can be estimated from the assumption of lift-to-drag ratio. Therefore, the design thrust in cruise at various Mach numbers can be estimated. For the Pratt & Whitney turbofan engines, thrust (T) is defined by EPR, so that the thrust model is set up as(3.23)T=f(h,W,M,CAS,EPR,ṁf).For GE turbofan engines, the rpm of the low-pressure compressor (N1) is used to set the level of thrust, so that the thrust model is set up as(3.24)T=f(h,W,M,CAS,N1,ṁf).In the present study, the P & W turbofan engines powering one twin-jet transport will be illustrated. The actual thrust in operation is obtained by using the recorded variables in the FDR, in particular the fuel flow rates.In the present study, the P & W engines powering the twin-jet commercial aircraft will be illustrated.The following climb equation [13] is to be satisfied in the least-square sense over a 5-second internal:(3.25)WgdVdt=T-D-Wsinγ,(3.26)DW=DLcosγ.All these equations are still valid in descent with negative climb angles (γ). The above equations are further employed for parameter identification in the process of modeling.Once the thrust model is generated as a function ofh, W, M, CAS, EPR, and ṁf with the flight conditions of climbing, cruise, and descent, one can estimate the thrust magnitude by inserting these flight variables from the FDR into the model.


### 3.5. Fuzzy-Logic Unsteady 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. For longitudinal aerodynamics, the models are assumed to be of the form [6](3.27)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 (3.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 [8],(3.28)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

The twin-jet transport of the present study encountered clear-air turbulence in cruise flight at the altitude around 10,050 m. As a result, several passengers and cabin crews sustained injuries, because of which this event was classified as an accident. The present study was initiated to examine possible concepts of accident prevention in the future. The datasets used for the modeling are extracted from the FDR during turbulence encounter lasting for 92 seconds.The main aircraft geometric and inertial characteristics are taken to beW (take-off) = 1,431,800 N (321900 lb),S= 260 m2  (2798.7 ft2), c̅= 6.608 m (21.68 ft), and b= 44.827 m (147.08 ft),Ixx= 10,710,000 kg-m2 (7,899,900 slugs-ft2), and Iyy= 14,883,800 kg-m2 (10,978,000 slugs-ft2),Izz= 25,283,271 kg-m2 (18,648,470 slugs-ft2) and Ixz= 0.0 kg-m2.The required operational parameters in FDR dataset for generating aerodynamic data files are time (t), CAS, pressure altitude (h), roll attitude (ϕ), pitch attitude (θ), magnetic heading (ψ), longitudinal acceleration (ax), lateral acceleration (ay), vertical acceleration (az), angle of attack (α), aileron deflection (δa), elevator (δe), rudder (δr), stabilizer (δs), engine EPR, outside air temperature, wind speed, 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:(3.1)ϕ̇=p+qsinϕtanθ+rcosϕtanθ,(3.2)θ̇=qcosϕ-rsinϕ,(3.3)ψ̇=(qsinϕ+rcosϕ)secθ,(3.4)V̇=(ax-gsinθ)cosαcosβ+(ay+gsinϕcosθ)sinβ+(az+gcosϕcosθ)sinαcosβ,(3.5)α̇=[(az+gcosθcosϕ)cosα-(ax-gsinθ)sinα](Vcosβ)+q-tanβ(pcosα+rsinα),(3.6)β̇=cosβ(ay+gcosθsinϕ)V+psinα-rcosα-sinβ[(az+gcosθcosϕ)sinα-(ax-gsinθ)cosα]V,
where g is acceleration due to gravity, V is flight speed, β is sideslip angle, p is roll rate, q is pitch rate, and r is yaw rate in (3.1)–(3.6). Let the biases be denoted by bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,andbψ, respectively, for ax, ay, az, and so forth. These biases are estimated by minimizing the squared sum of the differences between the two sides of the above equations. These equations in vector form can be written as(3.7)z⇀̇=f⇀(x)=f⇀(xm-Δx),
where(3.8)z⇀=(V,α,β,θ,ϕ,ψ)T,(3.9)x⇀m=(ax,ay,az,p,q,r,V,α,β,θ,ϕ,ψ)T,(3.10)Δx⇀=(bax,bay,baz,bp,bq,br,bV,bα,bβ,bθ,bϕ,bφ)T,
where the subscript “m” indicates the measured or recorded values. The cost function is defined as(3.11)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, variables not present in the FDR, such as β, p, q, and r, are also estimated.The force and moment coefficients are obtained from the following flight dynamic equations [7, 10] about the airplane body axes:(3.12)max=Cxq̅S+Tx,(3.13)may=Cyq̅S+Ty,(3.14)maz=Czq̅S+Tz,(3.15)Clq̅Sb=Ixxṗ-Ixz(ṙ+pq)-(Iyy-Izz)qr,(3.16)Cmq̅Sc̅=Iyyq̇-Ixz(r2-p2)-(Izz-Ixx)rp-Tm,(3.17)Cnq̅Sb=Izzṙ-Ixz(ṗ-qr)-(Ixx-Iyy)pq,
where m is denoted as aircraft mass; q̅ is denoted as dynamic pressure; S is denoted as wing reference area; Cx, Cz, and Cm are denoted as longitudinal aerodynamic force and moment coefficients; Cy, Cl, and Cn are denoted as lateral-directional aerodynamic force and moment coefficients; Ixx, Iyy, and Izz are denoted as moments of inertia about x-, y-, and z-axes, respectively; Ixy, Ixz, and Iyz are denoted as products of inertia; Tx, Ty, Tz, and Tm are denoted as thrust terms about x-, y-, z-axes, and in equation of pitching moment, respectively, in (3.7)–(3.12).The above equations are used to determine all aerodynamic coefficients based on accelerometer readings (ax, ay, and az), Euler angles (ϕ, θ, and ψ), angular rates (p, q, and r), and thrusts (Tx, Ty, Tz, and Tm). The angular rates are estimated through compatibility analysis. Since thrusts were not measured during flight for most flight vehicles, those values and the effects on the forces and pitching moments in equations of (3.7), (3.8), (3.9), and (3.11) should be predicted by a thrust model.


## 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 is 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 following the classical unsteady aerodynamic theory of Theodorsen [11]. For lateral-directional aerodynamics, it is based on the time variation of roll angle [8].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 [7, 12] as follows:(3.18)α(t)=α̅+αcos(ωt+ϕ̅),(3.19)α̇(t)=-aωsin(ωt+ϕ̅),
where those terms on the left hand side of (3.13) and (3.14) 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 angular frequency (ω). These unknowns are calculated through an optimization method by minimizing the following cost function (least squares):(3.20)J=∑i=1n[αi-(α̅+acos(ωti+ϕ̅))]2+[α̇i-(α̅+aωsin(ωti+ϕ̅))]2.In (3.15), α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 (3.15). 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(3.21)k1=ωc̅V,
where c̅ is the mean chord length of wing airfoil section.The lateral-directional equivalent reduced frequency is defined as(3.22)k2=ωbV,
where b is the wing span.


## 3.4. Fuzzy-Logic Thrust Model

As shown before, the thrust terms appear in the force equations and the pitching moment equations (3.7)–(3.9) and (3.11). Since the values of thrust for aircraft in flight cannot be directly measured in the current state of the art, they are not recorded in the FDR. The manufacturers of airplanes and engines and safety agencies agreed that using such parameters as the Mach number, airspeed, flight altitude, temperature, the rpm of the pressure compressors, and engine pressure ratios is adequate to estimate the engine thrust. 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 realistic one tied to the recorded engine performance parameters will be developed with the fuzzy-logic algorithm.For a commercial aircraft, most likely only the axial force and the pitching moment are affected by thrust. This assumption will be made in this paper. Theoretically, clear-air turbulence (i.e., random change inu, w (or α), and v (or β) affects the engine performance through its effects on static and dynamic distortions at the engine face. However, the effects of clear-air turbulence on the engine parameters and the thrust magnitude are not known and, therefore, are ignored.For this purpose, data from the flight manual for the fuel flow rates (ṁf) at various altitudes (h), weights (W), Mach numbers (M), calibrated airspeed (CAS), engine pressure ratios (EPRs) in cruise flight are utilized. Note that the drag polar for a given aircraft is generally not known to most researchers. To estimate it and, hence, the thrust magnitude in cruise, the assumption of a design lift-to-drag ratio (L/D) of 17.5 is made. This value of lift-to-drag in cruise is assumed based on the past design experience for twin-jet transports. In the flight manual, various weights, altitudes, Mach numbers, CAS, EPR, and fuel flow rates in cruise are tabulated. The lift coefficient can be calculated at each flight condition immediately. As a result, the drag coefficient can be estimated from the assumption of lift-to-drag ratio. Therefore, the design thrust in cruise at various Mach numbers can be estimated. For the Pratt & Whitney turbofan engines, thrust (T) is defined by EPR, so that the thrust model is set up as(3.23)T=f(h,W,M,CAS,EPR,ṁf).For GE turbofan engines, the rpm of the low-pressure compressor (N1) is used to set the level of thrust, so that the thrust model is set up as(3.24)T=f(h,W,M,CAS,N1,ṁf).In the present study, the P & W turbofan engines powering one twin-jet transport will be illustrated. The actual thrust in operation is obtained by using the recorded variables in the FDR, in particular the fuel flow rates.In the present study, the P & W engines powering the twin-jet commercial aircraft will be illustrated.The following climb equation [13] is to be satisfied in the least-square sense over a 5-second internal:(3.25)WgdVdt=T-D-Wsinγ,(3.26)DW=DLcosγ.All these equations are still valid in descent with negative climb angles (γ). The above equations are further employed for parameter identification in the process of modeling.Once the thrust model is generated as a function ofh, W, M, CAS, EPR, and ṁf with the flight conditions of climbing, cruise, and descent, one can estimate the thrust magnitude by inserting these flight variables from the FDR into the model.


## 3.5. Fuzzy-Logic Unsteady 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. For longitudinal aerodynamics, the models are assumed to be of the form [6](3.27)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 (3.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 [8],(3.28)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

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 (SSEs) and the square of multiple correlation coefficients (R2). Figure 2 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 have good agreement with the flight data. The Cm-data scattering is most likely caused by turbulence-induced buffeting on the structure, in particular on the horizontal tail. Once the aerodynamic models are set up, one can calculate all necessary derivatives to analyze the stability.Predicted aerodynamic coefficients in normal force and moments for a twin-jet transport encountering severe clear-air turbulence at cruise altitudes around 10,050 m.
(a)(b)(c)(d)The fuzzy-logic aerodynamic models are capable of generating the continuous derivatives for the static and dynamic stability study of a twin-jet transport in turbulence response. Firstly, how the fuzzy-logic prediction is achieved will be illustrated with one numerical example in theCz calculation. At first, the range for each variable is defined to be larger than what actually occurred in the present set of Cz  data as follows: [α]=[-13,12], [α̇]=[54,50], [q]=[-20,10], [k1]=[0,0.6], [β]=[-7,3], [δe]=[-10,6], [M]=[0,1.6], [p]=[-24,38], [δs]=[-3,3], and [q̅]=[4.964,21.746].For the first cell(1,1,1,1,1,1,1,1,1,1), the coefficients in (2.1) after model training are found to be pk1= (2.61755, 1.26662, 1.42338, 2.07962, -0.44241, 2.78017, 1.78150, 1.30818, 1.82872, 1.67592, 1.13787).Assume that in the following flight conditionsCz is to be predicted: α=6.91015 deg.; α̇=2.95510 deg/sec; q=1.16609 deg/sec; k1=0.01965; β=-1.55252; δe=0.68120 deg; M=0.77279; p=-2.62359 deg/sec; δs=-0.13930 deg; q̅=11.0545 kpa.These values of variables are converted to[0,1]. For example, xα=[6.91015-(-13)]/[12-(-13)]=0.79641.Other variable values are converted in the same way. It follows that the cell internal function becomesP1=2.61755+(1.26662)*(0.79641)+(1.42338)*(0.54764)+(2.07962)*(0.70554)-(0.44241)*(0.03275)+(2.78017)*(0.54475)+(1.7815)*(0.66758)+(1.30818)*(0.48299)+(1.82872)*(0.34478+(1.67592)*(0.47678)+(1.13787)*(0.3730)=11.04817.The membership functions for the first cell are exactly equal toxr, being 0.79641, 0.54764, and so forth. Their product can be calculated to be 1.08536E- 004. Therefore, the contribution of the first cell to the total output is 11.04817 * 1.08536E-004= 1.19912E-003.The total output from all cells can be calculated to be 5.9962; while the denominator in (2.6) is calculated to be 7.46459. Therefore, the final prediction is 0.8033. Comparing with data of 0.81038, this prediction has an error of -0.88%.To examine the stability characteristics, it is imperative to understand the flight environment in detail. The corresponding flight data are presented in Figure3. In general, one uses an instead of az for turbulence intensity study. Note an is the same as az in this paper. The variation of normal acceleration is presented in Figure 3(a), showing the highest an being 1.75 g around t=3930 seconds and the lowest being 0.02 g around t=3932 seconds. Figure 3(b) shows that α is approximately in phase with an. When an is the highest (around t=3930 seconds), the aircraft rapidly plunging downward with the altitude (h) reaching the lowest as shown in Figure 3(c); and α is highest about 6.5 deg. in Figure 3(b). At the same time, M is around 0.77 in Figure 3(d). Since α reaches a value about 6.5 deg in 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 is included in the total α.The time history of flight variables for a twin-jet transport in sever clear-air turbulence at the altitude around 10,050 m in transonic flight.
(a)(b)(c)(d)The aerodynamic derivatives extracted from the unsteady aerodynamic models can be calculated with two approaches, first one being the central difference method [7] and the other being the forced small-amplitude oscillation method [14]. The longitudinal stability derivative (Cmα) is extracted from the model of Cm. It is evaluated with the central difference approach as follows:(4.1)Cmα=[Cm(α+Δα,…)-Cm(α-Δα,…)]2Δα,
where Δα=0.5 degree represents that αis perturbed by 0.5 degree while keeping all other variables unchanged. The roll damping (Clp) is extracted from the models of Cl with the central difference approach as follows:(4.2)Clp=[Cl(…,p+Δp,…)-Cl(…,p-Δp,…)]2Δp,
where Δp deg/sec. Similarly, all other aerodynamic derivatives are calculated by using the same method.The time period between 3927.5 seconds and 3932.5 seconds is emphasized in evaluating the stability characteristics, because of the plunging motion that affects the flight safety the most. In order to evaluate the variations in characteristics, the units of all aerodynamic derivatives are converted to rad-1. The main longitudinal and lateral-directional stability derivatives along the flight path are presented in Figure 4. It should be noted that these derivatives are evaluated at the instantaneous conditions, instead of about the trim conditions as have been traditionally done. From the point of view in static stability, initially, the configuration has longitudinal stability (Czα<0 and Cmα<0) as shown in Figure 4(a), stable longitudinal damping (Cmq<0) in Figure 4(b), lateral stability (Clβ<0) and directional stability (Cnβ>0) in Figure 4(c), small roll damping (Clp<0) and insufficient directional damping (Cnr small or positive) in Figure 4(d). During the plunging motion, in the period between t = 3928.5 seconds. and t = 3930.5 seconds, Cmα>0 and Clβ>0, so that the static stability becomes unstable. 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 of the static stability derivatives along the flight path.
(a)(b)(c)(d)Figure5 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 5(a), the oscillatory derivatives are defined asThe time history of main longitudinal and lateral-directional oscillatory derivatives along the flight path.
(a)(b)(c)(d)(4.3)(Cmq)osc=Cmq+Cmα̇,(4.4)(Czq)osc=Czq+Czα̇In Figure4(c), the oscillatory derivatives are defined as(4.5)(Clp)osc=Clp+Clβ̇sinα,(4.6)(Cnr)osc=Cnr-Cnβ̇cosα.During the plunging motion, the values have some differences between oscillatory and damping derivatives in Figures5(a) and 5(c) due to the effects of the dynamic derivatives (i.e., α̇ and β̇-derivatives). 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 after t = 3929.5 seconds; while the effects of β̇-derivative is to cause the directional stability more unstable (i.e., (Cnr)osc more positive). These results indicate that the turbulent crosswind has the effects on directional stability and damping. Although the dynamic derivatives tend to be small for the present configuration, these are much helpful to understand the unknown factors of instability characteristics. 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.


## 5. Concluding Remarks

The main objective in this paper was to illustrate the nonlinear unsteady aerodynamic models based on the FLM technique having the capability to evaluate the variations in stability of commercial aircraft with adverse weather effects. It was shown that the FLM technique was capable of handling nonlinear and unsteady aerodynamic environment exhibited for a twin-jet transport in severe clear-air turbulence with sudden plunging motion in transonic flight. The predicted results showed that the models could produce relatively accurate aerodynamic coefficients and several derivatives for the assessment of stability characteristics, especially for the sensitive study of unknown factors in adverse weather conditions.



---

*Source: 924275-2010-03-10.xml*