# Structural Health Monitoring of Transport Aircraft with Fuzzy Logic Modeling

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

---

## Abstract

A structural health monitoring method based on the concept of static aeroelasticity is presented in this paper. This paper focuses on the estimation of these aeroelastic effects on older transport aircraft, in particular the structural components that are most affected, in severe atmospheric turbulence. Because the structural flexibility properties are mostly unknown to aircraft operators, only the trend, not the magnitude, of these effects is estimated. For this purpose, one useful concept in static aeroelastic effects for conventional aircraft structures is that under aeroelastic deformation the aerodynamic center should move aft. This concept is applied in the present paper by using the fuzzy-logic aerodynamic models. A twin-jet transport aircraft in severe atmospheric turbulence involving plunging motion is examined. It is found that the pitching moment derivatives in cruise with moderate to severe turbulence in transonic flight indicate some degree of abnormality in the stabilizer (i.e., the horizontal tail). Therefore, the horizontal tail is the most severely affected structural component of the aircraft probably caused by vibration under the dynamic loads induced by turbulence.


---

## Body

## 1. Introduction

The transport aircraft experiences frequently high loads, especially in severe turbulence encounters, sudden evasive motion in avoiding air collision, and hard landings. These high loads may cause structural damage, fatigue cracking, or nontypical structural deformation. These structural problems are most noticeable for aging transport aircraft. Structural integrity of an aging aircraft is determined by the level of fatigue, corrosion, and corrosion-assisted fatigue [1]. Each certificate holder of aircraft maintenance factory has the maintenance program for the aircraft that includes damage tolerance-based inspections and procedures. In conventional maintenance program, there are no perfect and sure ways to detect or predict the structural integrity of aging aircraft.In general, the airlines use nondestructive inspection (NDI) technology to determine structural faults of transport aircraft. However, NDI is very time consuming and expensive. On the other hand, if a structural analysis tool is to be employed, structural flexibility characteristics for an aging aircraft would be needed; but are usually not available. In addition, the aircraft structures may corrode in a faster rate than expected, and structural material degradation problems may also occur, especially when the aircraft is loaded near or beyond its design load limits [2]. As a result, the structural flexibility characteristics may not be the same as in the original design after a certain period of services.Dynamic aeroelasticity involves the analysis of structural response to dynamic external aerodynamic loads and aircraft flutter onset [3]. In the present paper, the focus is on static aeroelasticity that estimates the effects of aerodynamic loads on structural deformation which in turn changes the aerodynamic loads [4], as well as structural response, to dynamic loads.Cases of structure failure in midair are rare but are possible. In 2002, a transport aircraft disintegrated into four pieces while climbing to cruise altitude at FL350 (10,667 m). According to the final investigation report on the cause of the crash, the aircraft suffered a tail strike in a previous incident twenty-two years before [5]. This flight accident was the result of “metal fatigue” due to inadequate maintenance. Based on the same analysis to be performed in this paper, the empennage of this aircraft has had a serious structural deformation under high dynamic pressures since the previous flight [6]. Any empennage with fatigue deformation would change the effective setting of pitch trim and pitch controls. That is, the stabilizer and/or elevator would have abnormal offset, and flight control power has unusual changes. In other words, if the airlines could detect such structural fatigue and fix the problem, this mishap could have been avoided.A method based on the concept of aeroelastic effects on the aerodynamic center location and pitch trimming power to indirectly determine the structural integrity of transport aircraft is presented. This method is based on the flight data of the flight data recorder (FDR) in an accident flight or the quick access recorder (QAR) in a nominal flight. Aerodynamic models are established through fuzzy-logic modeling (FLM) of the flight data. In this paper, a twin-jet transport is employed for this study. Based on the service age, this aircraft is classified as an aging aircraft. This twin-jet transport encountered severe atmospheric turbulence with sudden plunging motion in a revenue flight.


## 2. Numerical Method Development

Since the aerodynamic models are established by using flight data, modeling technique is important and need to be carefully considered. First of all, these data are filtered through compatibility analysis [7] to satisfy the kinematic equations to remove measurement noises and biases, and aerodynamic coefficients are calculated from the flight dynamic equations. Modeling procedures start from setting up numerical relations between the input (i.e., flight dynamic variables) and output (i.e., aerodynamic coefficients). In order to obtain continuous derivatives, the present paper uses internal functions, instead of fuzzy sets, to generate the output of the model [8].System identification in the present paper includes two tasks: one is the model structure identification, and the other one is to identify the parameters that represent their corresponding model structures. The present modeling method was first developed by Takagi and Sugeno [9] and Sugeno and Tanaka [10], and later in 1995, Xie et al. [11, 12] applied the theory to simulate microelectronic processes with very good accuracy. The present paper is based on the modeling technique suggested by Xie et al. This technique was first applied to aviation technology in 1997 and later to flight data [13–16]. A new application of this modeling technique to examine aeroelastic effects of a transport aircraft is presented in this paper. A detailed description of the fuzzy-logic algorithm is available in the work of Chang et al. [7] and is summarized in Section 2.1.


### 2.1. Fuzzy-Logic Technique

The FLM method takes advantage of correlating multiple parameters without assuming explicit functional relations among them. The algorithm employs many internal functions to represent the contributions of fuzzy cells (to be defined later) to the overall prediction. The internal functions are assumed to be linear functions of input parameters as follows [7]:
(1)Pi=yi(x1,x2,…,xr,…,xk)=p0i+p1ix1+⋯+prixr+⋯+pkixk,
wherepri, r=0,1,2,…,k, are the coefficients of internal functions yi and k is number of input variables. In (1), yi is the 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.With the internal function chosen in a linear form, the fuzzy-logic model resembles the multiple linear regression method. What makes the fuzzy-logic model unique is that it is in a form of fuzzy cell structure composed of linear equations. In other words, there are different numbers of cells corresponding with each input parameter. The values of each input variable in the modeling, such as the angle of attack, are divided into several ranges, each of which representing a membership function withA(xr) as its membership grade. 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 [7] is stated as follows: if x1 isA1i(x1), x2 isA2i(x2),…, andxkisAki(xk), then the cell output can be stated as (1), 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 parametersp0i,p1i,…,pri,…,pki to be determined, andAki(xk)denotes the membership function for xk. Each function covers a certain range of input variables.For a given system with input variablesx1,x2,…,xr,…,xk of one data point, the recorded values of each input variables are normalized by using
(2)xr,norm=xr-xr,minxr,max-xr,min,r=1,2,…,k.Hereafterxr,normis denoted byxrfor simplicity in description. The membership grading ranges are transformed into the domain [0,1]: “0” meaning no effect from the corresponding internal function and “1” meaning full effect. Generally, overlapped triangles are frequently the shapes used to represent the grades. A fuzzy cell is formed by taking one membership function from each variable. The output of the fuzzy-logic model is the weighted average of all cell outputs.In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function, (1). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm. The output is estimated by the center of gravity method. For the jth input (x1,j,x2,j,…,xr,j,…,xk,j), the output is as follows:
(3)y^=∑i=1nop[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]Pi∑i=1nop[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]ddddddddddddd000d11ddddddddddj=1,2,…,m.
In (3), op[A1i(x1,j),…,Aki(xk,j)] is the weighted factor 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 records. The symbol “op” stands for product operator of its elements in the present paper.There are two main tasks involved in the fuzzy-logic modeling process. One is the determination of coefficients of the linear internal functions. The other is to identify the best structure of fuzzy cells of the model, that is, to determine the best number of membership functions for each input variable in the modeling. The coefficients are calculated with the gradient-descent method by minimizing the sum of squared errors (SSEs) [14]:
(4)SSE=∑j=1m(y^j-yj)2.On the other hand, the structure of fuzzy cells is optimized by maximizing the multiple correlation coefficients (R2):
(5)R2=1-{∑j=1m(y^j-yj)2}{∑j=1m(y--yj)2},
where y^j is the output of the fuzzy-logic model, yjis the measured data, and y- is the average value of all data.The aerodynamic model is defined by the values ofpri-coefficients. These coefficients are determined by minimizing SSE (4) with respect to these coefficients. Minimization is achieved by the gradient-descent method with an iterative formula defined by
(6)pr,t+1i=pr,ti-αr∂(SSE)∂pri,
where αr is convergence factor or step size in the gradient method; subscript index t denotes the iterative sequence.After simplification, (4) becomes  (7a)∂(SSE)∂pri=2∑j=1m(y^j-yj)∂y^j(x1,j,…,xk,j,pr1,…,pkn)∂pri.Since the computed gradient tends to be small with (7a) and the convergence is slow, to accelerate the convergence, the iterative formulas are modified by using the local squared errors to give the following.Forr=0,
(7b)p0,t+1i=p0,ti-2α0(y^j-yj)=×op[A1i(x1,j),…,Aki(xk,j)]∑s=1nop[A1s(x1,j),…,Aks(xk,j)]
and for r=1,…,k,
(7c)pr,t+1i=pr,ti-2αr(y^j-yj)=×op[A1i(x1,j),…,Aki(xk,j)]xr,j∑s=1nop[A1s(x1,j),…,Aks(xk,j)].The iteration during the search sequence stops when one of the following three criteria [7] is satisfied:(1)
cost =SSEt<ε1,(2)
RER =(SSEt-SSEt-1)/SSEt<ε2,(3)
t=tmax.In the above criteria,SSEI is the sum of squared errors (SSEs) being the cost function, and RER = (cost_current − cost_previous)/cost_current is 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. Figure 1 is the flowchart of parameter identification algorithm.Figure 1
Flowchart of parameter identification algorithm.Given membership functions and the training data, this parameter identification procedure can be applied to establish a fuzzy-logic model. Once the aerodynamic models are set up, one can input influencing variables and flight conditions to estimate aerodynamic coefficients. By varying the influencing parameters one at a time, the corresponding changes in the aerodynamic coefficients are used to predict the related stability and control derivatives for the analysis of flight dynamic characteristics. All aerodynamic derivatives are computed with central differences [15, 16].


### 2.2. Static Aeroelasticity

Since nondestructive inspection is too time-consuming and expensive to conduct frequently, it would be more economical and effective if a complementary computer method based on flight data in the digital FDR can be used. This method is based on the concept of static aeroelasticity as described in [17]. The idea is that any structural deformation in flight due to flight loads will change some aerodynamic parameters with a certain trend under normal conditions. The latter includes applicability of beam theory to deformation of all structural components, and the variation of bending and torsional rigidity can be represented by the average values of transport aircraft. Note that the beam theory is applied to the elastic axis of the aircraft structure and should be applicable to most transport aircrafts with metal structure. The aerodynamics was represented by an aerodynamic panel method in [17]. The formulation is summarized as follows. Note that although the following discussion is based on a panel or finite-element method, the results to be presented will be obtained by modeling the FDR data, not the panel method. The latter is only employed to show the concept.If an aircraft surface is divided intoN panels, then the panel forces would be
(8){FA}=q-[A]{αw+αd},
where q- is dynamic pressure; [A] is the aerodynamic influence coefficient in matrix form defined as the aerodynamic load per unit dynamic pressure per unit area at one panel due to a unit angle of attack at some other panels, {αd} is angle-of-attack change due to structural deformation in vector matrix form; {αw} is wing angle of attack in vector matrix form. The aerodynamic forces acting on the structure will cause the structural deformation. The angle-of-attack changes due to structural deformation in vector matrix form is as follows:
(9){αd}=[Cθ]{FA},
where Cθ is the so-called structural flexibility matrix, being defined as the local angle-of-attack change at one panel due to a unit load at some other panels.In general, the wing angle of attack is defined to be constant along the chordwise direction, so the vector matrix form can be simplified as(10){αw}={1}α,
where α is the magnitude of the aircraft angle of attack. Substituting (9) into (8) and solving for FA then give
(11){FA}=[I-q-[A][Cθ]]-1q-[A]{1}α,
where [I] is the identity matrix. Summing overall panels, the total force coefficient can be obtained as follows:
(12)Czαe={FA}{S}Tq-Sref=[I-q-[A][Cθ]]-1Czα,
where Czα=[A]{1}{S}T/Sref and S is a vector consisting of panel areas.Therefore, the elastic derivative(13)∂2Cze∂α∂q-=[I-q-[A][Cθ]]-1[A][Cθ][I-q-[A][Cθ]]-1Czα
can be used to demonstrate the magnitude of aeroelastic effects, mainly of the wing and tail. Note that the magnitude is proportional to the structural flexibility matrix [Cθ] which is usually very small in magnitude.Similarly, the derivative(14)∂2Cme∂α∂q-
indicates the α-slope change in pitching moment coefficient due to dynamic pressure and hence the structural deformation and may indicate the relative magnitude of empennage aeroelastic effects. The reason for using the second derivatives is that the first derivatives with dynamic pressure are not sensitive enough to exhibit the changes and tend to be too small. These derivatives are all calculated with a central difference formula by using the numerical models. Note that the aerodynamic center (ac) location is determined by the following relation:
(15)x-ac-x-cg=-∂Cm∂CL,
where x-is the length measured from the leading edge of the mean aerodynamic chord to the aerodynamic center (ac) or center of gravity (cg) and is nondimensionalized with the mean aerodynamic chord. Therefore, if the second derivative of Cme is negative, x-ac would be increased, and hence the aerodynamic center is moved back in accordance with the estimation of aeroelastic effects in reference [17]. In addition, the wing should twist nose down, so that the Cze derivative could be negative for a normal wing structure.


## 2.1. Fuzzy-Logic Technique

The FLM method takes advantage of correlating multiple parameters without assuming explicit functional relations among them. The algorithm employs many internal functions to represent the contributions of fuzzy cells (to be defined later) to the overall prediction. The internal functions are assumed to be linear functions of input parameters as follows [7]:
(1)Pi=yi(x1,x2,…,xr,…,xk)=p0i+p1ix1+⋯+prixr+⋯+pkixk,
wherepri, r=0,1,2,…,k, are the coefficients of internal functions yi and k is number of input variables. In (1), yi is the 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.With the internal function chosen in a linear form, the fuzzy-logic model resembles the multiple linear regression method. What makes the fuzzy-logic model unique is that it is in a form of fuzzy cell structure composed of linear equations. In other words, there are different numbers of cells corresponding with each input parameter. The values of each input variable in the modeling, such as the angle of attack, are divided into several ranges, each of which representing a membership function withA(xr) as its membership grade. 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 [7] is stated as follows: if x1 isA1i(x1), x2 isA2i(x2),…, andxkisAki(xk), then the cell output can be stated as (1), 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 parametersp0i,p1i,…,pri,…,pki to be determined, andAki(xk)denotes the membership function for xk. Each function covers a certain range of input variables.For a given system with input variablesx1,x2,…,xr,…,xk of one data point, the recorded values of each input variables are normalized by using
(2)xr,norm=xr-xr,minxr,max-xr,min,r=1,2,…,k.Hereafterxr,normis denoted byxrfor simplicity in description. The membership grading ranges are transformed into the domain [0,1]: “0” meaning no effect from the corresponding internal function and “1” meaning full effect. Generally, overlapped triangles are frequently the shapes used to represent the grades. A fuzzy cell is formed by taking one membership function from each variable. The output of the fuzzy-logic model is the weighted average of all cell outputs.In each fuzzy cell, the contribution to the outcome (i.e., the cell output) is based on the internal function, (1). The final prediction of the outcome is the weighted average of all cell outputs after the process of reasoning algorithm. The output is estimated by the center of gravity method. For the jth input (x1,j,x2,j,…,xr,j,…,xk,j), the output is as follows:
(3)y^=∑i=1nop[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]Pi∑i=1nop[Ai(x1,j),…,Ai(xr,j),…,Ai(xk,j)]ddddddddddddd000d11ddddddddddj=1,2,…,m.
In (3), op[A1i(x1,j),…,Aki(xk,j)] is the weighted factor 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 records. The symbol “op” stands for product operator of its elements in the present paper.There are two main tasks involved in the fuzzy-logic modeling process. One is the determination of coefficients of the linear internal functions. The other is to identify the best structure of fuzzy cells of the model, that is, to determine the best number of membership functions for each input variable in the modeling. The coefficients are calculated with the gradient-descent method by minimizing the sum of squared errors (SSEs) [14]:
(4)SSE=∑j=1m(y^j-yj)2.On the other hand, the structure of fuzzy cells is optimized by maximizing the multiple correlation coefficients (R2):
(5)R2=1-{∑j=1m(y^j-yj)2}{∑j=1m(y--yj)2},
where y^j is the output of the fuzzy-logic model, yjis the measured data, and y- is the average value of all data.The aerodynamic model is defined by the values ofpri-coefficients. These coefficients are determined by minimizing SSE (4) with respect to these coefficients. Minimization is achieved by the gradient-descent method with an iterative formula defined by
(6)pr,t+1i=pr,ti-αr∂(SSE)∂pri,
where αr is convergence factor or step size in the gradient method; subscript index t denotes the iterative sequence.After simplification, (4) becomes  (7a)∂(SSE)∂pri=2∑j=1m(y^j-yj)∂y^j(x1,j,…,xk,j,pr1,…,pkn)∂pri.Since the computed gradient tends to be small with (7a) and the convergence is slow, to accelerate the convergence, the iterative formulas are modified by using the local squared errors to give the following.Forr=0,
(7b)p0,t+1i=p0,ti-2α0(y^j-yj)=×op[A1i(x1,j),…,Aki(xk,j)]∑s=1nop[A1s(x1,j),…,Aks(xk,j)]
and for r=1,…,k,
(7c)pr,t+1i=pr,ti-2αr(y^j-yj)=×op[A1i(x1,j),…,Aki(xk,j)]xr,j∑s=1nop[A1s(x1,j),…,Aks(xk,j)].The iteration during the search sequence stops when one of the following three criteria [7] is satisfied:(1)
cost =SSEt<ε1,(2)
RER =(SSEt-SSEt-1)/SSEt<ε2,(3)
t=tmax.In the above criteria,SSEI is the sum of squared errors (SSEs) being the cost function, and RER = (cost_current − cost_previous)/cost_current is 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. Figure 1 is the flowchart of parameter identification algorithm.Figure 1
Flowchart of parameter identification algorithm.Given membership functions and the training data, this parameter identification procedure can be applied to establish a fuzzy-logic model. Once the aerodynamic models are set up, one can input influencing variables and flight conditions to estimate aerodynamic coefficients. By varying the influencing parameters one at a time, the corresponding changes in the aerodynamic coefficients are used to predict the related stability and control derivatives for the analysis of flight dynamic characteristics. All aerodynamic derivatives are computed with central differences [15, 16].


## 2.2. Static Aeroelasticity

Since nondestructive inspection is too time-consuming and expensive to conduct frequently, it would be more economical and effective if a complementary computer method based on flight data in the digital FDR can be used. This method is based on the concept of static aeroelasticity as described in [17]. The idea is that any structural deformation in flight due to flight loads will change some aerodynamic parameters with a certain trend under normal conditions. The latter includes applicability of beam theory to deformation of all structural components, and the variation of bending and torsional rigidity can be represented by the average values of transport aircraft. Note that the beam theory is applied to the elastic axis of the aircraft structure and should be applicable to most transport aircrafts with metal structure. The aerodynamics was represented by an aerodynamic panel method in [17]. The formulation is summarized as follows. Note that although the following discussion is based on a panel or finite-element method, the results to be presented will be obtained by modeling the FDR data, not the panel method. The latter is only employed to show the concept.If an aircraft surface is divided intoN panels, then the panel forces would be
(8){FA}=q-[A]{αw+αd},
where q- is dynamic pressure; [A] is the aerodynamic influence coefficient in matrix form defined as the aerodynamic load per unit dynamic pressure per unit area at one panel due to a unit angle of attack at some other panels, {αd} is angle-of-attack change due to structural deformation in vector matrix form; {αw} is wing angle of attack in vector matrix form. The aerodynamic forces acting on the structure will cause the structural deformation. The angle-of-attack changes due to structural deformation in vector matrix form is as follows:
(9){αd}=[Cθ]{FA},
where Cθ is the so-called structural flexibility matrix, being defined as the local angle-of-attack change at one panel due to a unit load at some other panels.In general, the wing angle of attack is defined to be constant along the chordwise direction, so the vector matrix form can be simplified as(10){αw}={1}α,
where α is the magnitude of the aircraft angle of attack. Substituting (9) into (8) and solving for FA then give
(11){FA}=[I-q-[A][Cθ]]-1q-[A]{1}α,
where [I] is the identity matrix. Summing overall panels, the total force coefficient can be obtained as follows:
(12)Czαe={FA}{S}Tq-Sref=[I-q-[A][Cθ]]-1Czα,
where Czα=[A]{1}{S}T/Sref and S is a vector consisting of panel areas.Therefore, the elastic derivative(13)∂2Cze∂α∂q-=[I-q-[A][Cθ]]-1[A][Cθ][I-q-[A][Cθ]]-1Czα
can be used to demonstrate the magnitude of aeroelastic effects, mainly of the wing and tail. Note that the magnitude is proportional to the structural flexibility matrix [Cθ] which is usually very small in magnitude.Similarly, the derivative(14)∂2Cme∂α∂q-
indicates the α-slope change in pitching moment coefficient due to dynamic pressure and hence the structural deformation and may indicate the relative magnitude of empennage aeroelastic effects. The reason for using the second derivatives is that the first derivatives with dynamic pressure are not sensitive enough to exhibit the changes and tend to be too small. These derivatives are all calculated with a central difference formula by using the numerical models. Note that the aerodynamic center (ac) location is determined by the following relation:
(15)x-ac-x-cg=-∂Cm∂CL,
where x-is the length measured from the leading edge of the mean aerodynamic chord to the aerodynamic center (ac) or center of gravity (cg) and is nondimensionalized with the mean aerodynamic chord. Therefore, if the second derivative of Cme is negative, x-ac would be increased, and hence the aerodynamic center is moved back in accordance with the estimation of aeroelastic effects in reference [17]. In addition, the wing should twist nose down, so that the Cze derivative could be negative for a normal wing structure.


## 3. Aircraft Aerodynamic Modeling

### 3.1. Flight Data for Modeling

Two datasets are employed for comparative analysis of static aeroelastic effects on longitudinal aerodynamics for the twin-jet transport before and in severe atmospheric turbulence. One dataset is the aircraft response data to turbulence recorded in the FDR. Figure2 presents the time history of main flight variables for the twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in cruise phase. The variations of flight altitudes and angles of attack are indicated in Figures 2(a) and 2(b), respectively. The detailed aerodynamic environment will be described in Section 4.1. The largest dynamic pressure is 11,972 N/m2 (250 lb/ft2) before encountering severe turbulence, and the magnitudes of dynamic pressure decreased during the period of encountering the severe turbulence in the period time of t=3927~3940 sec, as indicated in Figure 2(c).The time history of main flight variables for the twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in cruise phase.
(a)
(b)
(c)Usually, the datasets in the climb phase are more suitable for the study of static aeroelastic effects, because the aircraft passes through the high dynamic pressure region before reaching the cruise altitude. Figure3 presents the time history of main flight variables for the transport aircraft in the climb phase. The dynamic pressure is more than 14,367 N/m2 (300 lb/ft2) in the time period of t=280~860 sec, while the range of flight altitudes is 4,572 m~8,534 m.The time history of main flight variables for the transport aircraft in the climb phase.
(a)
(b)
(c)The main aircraft geometric and inertial characteristics are taken to beW (takeoff) = 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), Iyy=14,883,800 kg-m2 (10,978,000 slugs-ft2),Izz=25,283,271 kg-m2 (18,648,470 slugs-ft2), Ixz=0.0 kg-m2.These two datasets for modeling are extracted from the FDR of the same flight. The service age of this transport aircraft is 14 years, 10 months, and 14 days, which should be close to the service years of an aging aircraft. The necessary data through modeling to determine the aerodynamics are time(t), CAS, pressure altitude (h), roll attitude (ϕ), pitch attitude (θ), magnetic heading (ψ), normal acceleration (az), lateral acceleration (ay), longitudinal acceleration (ax), 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. Modeling Preparation

As indicated earlier, the measured data in FDR for generating aerodynamic data must be checked and adjusted to satisfy the kinematic equations before modeling. The measurement noise and biases can be accounted for through compatibility analysis [7]. The time derivatives are all calculated using central differences. As a result of the analysis, variables not present in the FDR, such as sideslip angle (β) and angular rates (p, q, and r), are estimated. In the regression method, the force and moment coefficients are obtained from the flight dynamic equations [7] along the airplane body axes.


### 3.3. Longitudinal Main 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 the aircraft response to structure deformations. The longitudinal main aerodynamics is assumed to depend on the following ten flight variables [14]:
(16)Cz,Cm=f(a,α˙,q,k1,β,δe,M,p,δs,q-),
where the left-hand side represents the coefficients of normal force (Cz) and pitching moment (Cm) respectively. The variables on the right hand side of (16), where q is the pitch rate, k1 the reduced frequency, δe the elevator angle, M the Mach number, p the roll rate, δs the stabilizer angle, and q- the dynamic pressure. Note that all derivatives are calculated with the established aerodynamic models at the instantaneous conditions, or the local slopes, not around the trim conditions as usually done in the conventional method. 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.


## 3.1. Flight Data for Modeling

Two datasets are employed for comparative analysis of static aeroelastic effects on longitudinal aerodynamics for the twin-jet transport before and in severe atmospheric turbulence. One dataset is the aircraft response data to turbulence recorded in the FDR. Figure2 presents the time history of main flight variables for the twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in cruise phase. The variations of flight altitudes and angles of attack are indicated in Figures 2(a) and 2(b), respectively. The detailed aerodynamic environment will be described in Section 4.1. The largest dynamic pressure is 11,972 N/m2 (250 lb/ft2) before encountering severe turbulence, and the magnitudes of dynamic pressure decreased during the period of encountering the severe turbulence in the period time of t=3927~3940 sec, as indicated in Figure 2(c).The time history of main flight variables for the twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in cruise phase.
(a)
(b)
(c)Usually, the datasets in the climb phase are more suitable for the study of static aeroelastic effects, because the aircraft passes through the high dynamic pressure region before reaching the cruise altitude. Figure3 presents the time history of main flight variables for the transport aircraft in the climb phase. The dynamic pressure is more than 14,367 N/m2 (300 lb/ft2) in the time period of t=280~860 sec, while the range of flight altitudes is 4,572 m~8,534 m.The time history of main flight variables for the transport aircraft in the climb phase.
(a)
(b)
(c)The main aircraft geometric and inertial characteristics are taken to beW (takeoff) = 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), Iyy=14,883,800 kg-m2 (10,978,000 slugs-ft2),Izz=25,283,271 kg-m2 (18,648,470 slugs-ft2), Ixz=0.0 kg-m2.These two datasets for modeling are extracted from the FDR of the same flight. The service age of this transport aircraft is 14 years, 10 months, and 14 days, which should be close to the service years of an aging aircraft. The necessary data through modeling to determine the aerodynamics are time(t), CAS, pressure altitude (h), roll attitude (ϕ), pitch attitude (θ), magnetic heading (ψ), normal acceleration (az), lateral acceleration (ay), longitudinal acceleration (ax), 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. Modeling Preparation

As indicated earlier, the measured data in FDR for generating aerodynamic data must be checked and adjusted to satisfy the kinematic equations before modeling. The measurement noise and biases can be accounted for through compatibility analysis [7]. The time derivatives are all calculated using central differences. As a result of the analysis, variables not present in the FDR, such as sideslip angle (β) and angular rates (p, q, and r), are estimated. In the regression method, the force and moment coefficients are obtained from the flight dynamic equations [7] along the airplane body axes.


## 3.3. Longitudinal Main 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 the aircraft response to structure deformations. The longitudinal main aerodynamics is assumed to depend on the following ten flight variables [14]:
(16)Cz,Cm=f(a,α˙,q,k1,β,δe,M,p,δs,q-),
where the left-hand side represents the coefficients of normal force (Cz) and pitching moment (Cm) respectively. The variables on the right hand side of (16), where q is the pitch rate, k1 the reduced frequency, δe the elevator angle, M the Mach number, p the roll rate, δs the stabilizer angle, and q- the dynamic pressure. Note that all derivatives are calculated with the established aerodynamic models at the instantaneous conditions, or the local slopes, not around the trim conditions as usually done in the conventional method. 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.


## 4. Numerical Results and Discussions

### 4.1. Aerodynamic Environment

To examine the aeroelastic effects on longitudinal aerodynamics for the twin-jet transport during severe atmospheric turbulence, it is imperative to understand the flight environment first. The corresponding flight data are presented in Figure4. The variation of normal acceleration (az) is presented in Figure 4(a), showing the largest az being 1.75 g around t=3930sec and the lowest being 0.02 g around t=3932sec. Figure 4(b) shows that α is approximately in phase with az, and α reaches about 7 deg, which is well above the cruise value in transonic flight. Therefore, 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 α. When the aircraft is rapidly plunging downward with the altitude (h) reaching the lowest as shown in Figure 4(c), az becomes the largest (around t=3930sec). The Mach number (M) drops from 0.81 to 0.72 in Figure 4(d). The aeroelastic effects can be expected to be very significant under the circumstances of high g loads and rapidly changing az in transonic flight.The time history of flight variables for a twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in transonic flight.
(a)
(b)
(c)
(d)


### 4.2. Analysis of Model Predictions

In the present study, the accuracy of the established unsteady aerodynamic models with two (may be six in another application) aerodynamic coefficients through FLM algorithm is estimated by the sum of squared errors (SSEs) and the square of multiple correlation coefficients (R2). The elastic derivatives in the study of structural integrity are calculated with these aerodynamic models of aerodynamic coefficients.In the climb phase, initial modeling shows that the squared correlation coefficient (R2) of the normal force Cz and pitching moment Cm models can only reach 0.9882 and 0.7862, respectively. This may imply that uncertainties due to structural vibration, the estimated thrust, or some other unknown factors exist. At the same time, the computing time can be minimized if the number of data points can be reduced. For this purpose, the data reduction by model-based filtering [13] is applied. Firstly, the model structure is determined with all data included. After R2 remaining unchanged, all data points that deviate from the model-predicted values by more than a fixed percentage, such as 10% or larger for the present Cz and Cm data, will be deleted. After a few more run, the fixed percentage is reduced, and so on, until R2 reaches a value greater than a preset value. According to the experience in aerodynamic modeling, the value of R2 for Cz is usually higher. The value in the climb phase reaches 0.99. The value of R2 for Cm is lower than that of Cz. To avoid too many points being removed, the target value of R2 is 0.95. Since the present modeling is based on the least-square method, the final estimated model will not be affected by this type of filtering. In the climb phase, R2 of Cz and Cm models reaches from 0.9882 to 0.9948 and from 0.7862 to 0.9593, respectively. For the cruise data in atmospheric turbulence, all data points with Δt=0.025 sec should be included to cover the effects of all possible frequencies in turbulence.The predicted normal forceCz and pitching moment Cm coefficients for the climb and cruise phases are presented in Figures 5 and 6, respectively. The squared correlation coefficients (R2) of Cz and Cm for these two phases are shown in Table 1.Table 1
Squared correlation coefficients of numerical models.
R
2 for model
Climb
Cruise
R
2 for Cz (%)
99.48
97.87
R
2 for Cm (%)
96.06
95.93Predicted longitudinal aerodynamic coefficients in climb phase based on FDR data.
(a)
(b)Predicted longitudinal aerodynamic coefficients in the cruise phase based on FDR data.
(a)
(b)Note thatR2=0.9948 implies that 99.48% of the data can be explained and represented by the established numerical model [18]. These high correlation coefficients are reflected in the good agreement of the model-predicted results with the data.


### 4.3. Analysis of Aeroelastic Effects

Figure7 presents the time history of aeroelastic effects on longitudinal aerodynamics in the climb phase within a time span of 900 sec. The aeroelastic effects on longitudinal aerodynamics in the cruise phase within the time span of 92 sec are presented in Figure 8.Aeroelastic effects on longitudinal aerodynamics in climb phase within the time span of 900 sec.
(a)
(b)
(c)Aeroelastic effects on longitudinal aerodynamics in cruise phase within the time span of 92 sec.
(a)
(b)
(c)Again, based on the concept in static aeroelastic effects, deformation under aerodynamic loads for a conventional normal aircraft structure would decreaseCzαq¯ and make Cmαq¯ more negative, in such a way that the aerodynamic center would shift back [17]. This is also related to aeroelastic stability. For example, if Czαq¯ is positive, it means the wing would twist up under load, and the load would further increase.In the climb phase, there were not much atmospheric turbulence and tail wind. As shown in Figure7, although the dynamic pressure is high (Figure 7(a)), the magnitude of Cmαq¯ is very small as indicated in Figure 7(c). However, there was some wing bending so that Czαq¯ is negative as shown in Figure 7(c). As indicated earlier, negative Czαq¯ is normal. On the other hand, in the cruise phase, there were strong tail wind and varying turbulence intensity. As shown in Figures 8(b) and 8(c), there is a large variation for both magnitudes of Czαq¯ and Cmαq¯ in the severe turbulence encountered in the time period of t= 3927~3940 sec. The magnitudes of both Czαq¯ and Cmαq¯ in the cruise phase (Figures 8(b) and 8(c)) are all higher than that in the climb phase (Figures 7(b) and 7(c)). The magnitude of Czαq¯ after the plunging motion becomes less negative, implying the wing structure having less twisted nose down. Nose-down twisting under high loadings is normal for a typical sweptback wing. On the other hand, since Cmαq¯ could be positive and negative, as shown in Figure 8(c), with a positive value being not normal, the not-normal contribution must be from the stabilizer, because the latter contributes the most to the pitching moment. This is verified in Figure 9, where the aeroelastic parameter, ∂2Cm/∂δs∂q-, generates some positive values similar to Figure 8(c). It may be concluded that more significant aeroelastic effects in cruise shown in Figures 8(b) and 8(c) are caused by the structural vibration which is, in turn, caused by turbulence-induced dynamic loads.Figure 9
Aeroelastic effects on stabilizer effectiveness derivative in cruise.It should be noted that a pilot is always advised not to try to control the airplane in severe turbulence to avoid inappropriate input. In the present case, autopilot was deactivated in the severe turbulence. Therefore, the deflection to −1 deg as shown in Figure10(a) must be caused by the resulting dynamic loads from turbulence. A deflection of −1 deg (nose down) by the stabilizer at a transonic speed is significant. It may be concluded that structural integrity of the stabilizer, including the jackscrew, should be inspected more carefully, in particular the end play of jackscrew. The jackscrew is used to raise or lower the leading edge of the stabilizer. As shown in Figure 10(b), there appears to be some free play of the screw thread as the deflection after the severe turbulence encounter did not recover to a positive value which existed before the encounter. Under high dynamic loads, the thread wear could accelerate. Failure of the thread has doomed one airplane in an operation before [19]. Measurement of the jackscrew end play is part of the requirements in maintenance to monitor the structural integrity.Stabilizer deflection in cruise in atmospheric turbulence.
(a)
(b)


## 4.1. Aerodynamic Environment

To examine the aeroelastic effects on longitudinal aerodynamics for the twin-jet transport during severe atmospheric turbulence, it is imperative to understand the flight environment first. The corresponding flight data are presented in Figure4. The variation of normal acceleration (az) is presented in Figure 4(a), showing the largest az being 1.75 g around t=3930sec and the lowest being 0.02 g around t=3932sec. Figure 4(b) shows that α is approximately in phase with az, and α reaches about 7 deg, which is well above the cruise value in transonic flight. Therefore, 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 α. When the aircraft is rapidly plunging downward with the altitude (h) reaching the lowest as shown in Figure 4(c), az becomes the largest (around t=3930sec). The Mach number (M) drops from 0.81 to 0.72 in Figure 4(d). The aeroelastic effects can be expected to be very significant under the circumstances of high g loads and rapidly changing az in transonic flight.The time history of flight variables for a twin-jet transport in severe atmospheric turbulence at the altitude around 10,050 m in transonic flight.
(a)
(b)
(c)
(d)


## 4.2. Analysis of Model Predictions

In the present study, the accuracy of the established unsteady aerodynamic models with two (may be six in another application) aerodynamic coefficients through FLM algorithm is estimated by the sum of squared errors (SSEs) and the square of multiple correlation coefficients (R2). The elastic derivatives in the study of structural integrity are calculated with these aerodynamic models of aerodynamic coefficients.In the climb phase, initial modeling shows that the squared correlation coefficient (R2) of the normal force Cz and pitching moment Cm models can only reach 0.9882 and 0.7862, respectively. This may imply that uncertainties due to structural vibration, the estimated thrust, or some other unknown factors exist. At the same time, the computing time can be minimized if the number of data points can be reduced. For this purpose, the data reduction by model-based filtering [13] is applied. Firstly, the model structure is determined with all data included. After R2 remaining unchanged, all data points that deviate from the model-predicted values by more than a fixed percentage, such as 10% or larger for the present Cz and Cm data, will be deleted. After a few more run, the fixed percentage is reduced, and so on, until R2 reaches a value greater than a preset value. According to the experience in aerodynamic modeling, the value of R2 for Cz is usually higher. The value in the climb phase reaches 0.99. The value of R2 for Cm is lower than that of Cz. To avoid too many points being removed, the target value of R2 is 0.95. Since the present modeling is based on the least-square method, the final estimated model will not be affected by this type of filtering. In the climb phase, R2 of Cz and Cm models reaches from 0.9882 to 0.9948 and from 0.7862 to 0.9593, respectively. For the cruise data in atmospheric turbulence, all data points with Δt=0.025 sec should be included to cover the effects of all possible frequencies in turbulence.The predicted normal forceCz and pitching moment Cm coefficients for the climb and cruise phases are presented in Figures 5 and 6, respectively. The squared correlation coefficients (R2) of Cz and Cm for these two phases are shown in Table 1.Table 1
Squared correlation coefficients of numerical models.
R
2 for model
Climb
Cruise
R
2 for Cz (%)
99.48
97.87
R
2 for Cm (%)
96.06
95.93Predicted longitudinal aerodynamic coefficients in climb phase based on FDR data.
(a)
(b)Predicted longitudinal aerodynamic coefficients in the cruise phase based on FDR data.
(a)
(b)Note thatR2=0.9948 implies that 99.48% of the data can be explained and represented by the established numerical model [18]. These high correlation coefficients are reflected in the good agreement of the model-predicted results with the data.


## 4.3. Analysis of Aeroelastic Effects

Figure7 presents the time history of aeroelastic effects on longitudinal aerodynamics in the climb phase within a time span of 900 sec. The aeroelastic effects on longitudinal aerodynamics in the cruise phase within the time span of 92 sec are presented in Figure 8.Aeroelastic effects on longitudinal aerodynamics in climb phase within the time span of 900 sec.
(a)
(b)
(c)Aeroelastic effects on longitudinal aerodynamics in cruise phase within the time span of 92 sec.
(a)
(b)
(c)Again, based on the concept in static aeroelastic effects, deformation under aerodynamic loads for a conventional normal aircraft structure would decreaseCzαq¯ and make Cmαq¯ more negative, in such a way that the aerodynamic center would shift back [17]. This is also related to aeroelastic stability. For example, if Czαq¯ is positive, it means the wing would twist up under load, and the load would further increase.In the climb phase, there were not much atmospheric turbulence and tail wind. As shown in Figure7, although the dynamic pressure is high (Figure 7(a)), the magnitude of Cmαq¯ is very small as indicated in Figure 7(c). However, there was some wing bending so that Czαq¯ is negative as shown in Figure 7(c). As indicated earlier, negative Czαq¯ is normal. On the other hand, in the cruise phase, there were strong tail wind and varying turbulence intensity. As shown in Figures 8(b) and 8(c), there is a large variation for both magnitudes of Czαq¯ and Cmαq¯ in the severe turbulence encountered in the time period of t= 3927~3940 sec. The magnitudes of both Czαq¯ and Cmαq¯ in the cruise phase (Figures 8(b) and 8(c)) are all higher than that in the climb phase (Figures 7(b) and 7(c)). The magnitude of Czαq¯ after the plunging motion becomes less negative, implying the wing structure having less twisted nose down. Nose-down twisting under high loadings is normal for a typical sweptback wing. On the other hand, since Cmαq¯ could be positive and negative, as shown in Figure 8(c), with a positive value being not normal, the not-normal contribution must be from the stabilizer, because the latter contributes the most to the pitching moment. This is verified in Figure 9, where the aeroelastic parameter, ∂2Cm/∂δs∂q-, generates some positive values similar to Figure 8(c). It may be concluded that more significant aeroelastic effects in cruise shown in Figures 8(b) and 8(c) are caused by the structural vibration which is, in turn, caused by turbulence-induced dynamic loads.Figure 9
Aeroelastic effects on stabilizer effectiveness derivative in cruise.It should be noted that a pilot is always advised not to try to control the airplane in severe turbulence to avoid inappropriate input. In the present case, autopilot was deactivated in the severe turbulence. Therefore, the deflection to −1 deg as shown in Figure10(a) must be caused by the resulting dynamic loads from turbulence. A deflection of −1 deg (nose down) by the stabilizer at a transonic speed is significant. It may be concluded that structural integrity of the stabilizer, including the jackscrew, should be inspected more carefully, in particular the end play of jackscrew. The jackscrew is used to raise or lower the leading edge of the stabilizer. As shown in Figure 10(b), there appears to be some free play of the screw thread as the deflection after the severe turbulence encounter did not recover to a positive value which existed before the encounter. Under high dynamic loads, the thread wear could accelerate. Failure of the thread has doomed one airplane in an operation before [19]. Measurement of the jackscrew end play is part of the requirements in maintenance to monitor the structural integrity.Stabilizer deflection in cruise in atmospheric turbulence.
(a)
(b)


## 5. Concluding Remarks

The main objective of this paper was to present an evaluating method of aeroelastic effects without having the values of structural flexibility matrix to solve the mathematical problems in engineering. The method was based on the flight data extracted from the FDR to predict the aeroelastic effects for a twin-jet transport in severe atmospheric turbulence. The aeroelastic effects of the wing and empennage in severe atmospheric turbulence were significant due to the structural deformation under high dynamic loads in transonic flight and were found to be larger than that in the climb phase without turbulence effects. Because this method employed the flight data stored in the FDR, nondestructive inspection could be avoided in cases with negligible aeroelastic effects. In essence, the normal force and pitching moment of second derivatives were used to estimate qualitatively the aeroelastic effects on the wing and the tail, respectively. This method could be used to assist airlines to monitor the structural integrity as a complementary tool to improve aviation safety and operational efficiency.



---

*Source: 640852-2013-06-12.xml*