The E ﬀ ects of Icing on Aircraft Longitudinal Aerodynamic Characteristics

: To study the e ﬀ ects of ice accretion on the longitudinal aerodynamic characteristics of an aircraft, a two-part method for predicting longitudinal aerodynamic derivatives of iced aircraft is proposed. For the aircraft with a ﬂight test, a parameter identiﬁcation system based on maximum likelihood criterion and a longitudinal nonlinear ﬂight dynamics model is established. For the aircraft without a ﬂight test, an engineering prediction method of aerodynamic derivatives based on an individual component CFD calculation and narrow strip theory is established. According to the ﬂight test data of DHC-6 Twin Otter aircraft from NASA, the longitudinal aerodynamic parameters of both clean and artiﬁcially iced aircraft are obtained. Additionally, the longitudinal aerodynamic derivatives of the iced aircraft are calculated. Then, the correctness of the prediction method is veriﬁed by comparing the calculated results with the identiﬁcation results. The comparison of these results shows that the prediction method is correct and accurate, and it can be used to calculate the e ﬀ ects of icing on the aircraft longitudinal aerodynamic parameters.


Introduction
Aircraft icing is the phenomena of ice accretion on the aircraft. Aircraft icing will change the air flow around the lift surface, thus reducing the performance and control ability of the aircraft. Severe icing will cause air separation on the airfoil ahead of time, resulting in a stall [1], which poses a threat to the safety of flight characteristics. In order to prevent the serious consequences of icing, an anti-icing system is usually installed on the important lifting components of modern aircraft. However, flights with ice accretion on the aircraft are still unavoidable, and the accidents are usually caused by the failure or improper operation of ice protection systems. According to the data given by National Transportation Safety Board (NTSB) [2], from 1981 to 1988, about 542 aircraft accidents are caused by aircraft icing. According to the statistics of the International Civil Aviation Organization (ICAO), the death rate of passengers in aircraft icing accidents reached 39% in the early 1990s [3]. Aircraft icing is still a big problem for flight safety and the studies on this area are attracting more attention.
In order to define the effects of ice accretion on aircraft, the methods of parameter identification were used to extract the aerodynamic parameters from flight test data. Ranaudo et al. [4] compared the stability and control derivatives between clean and naturally iced aircraft, using a modified maximum likelihood estimation method with flight test data of Twin Otter. Ratvasky and Ranaudo [5] studied the effects of ice accretion on Twin Otter stability and control derivatives with flight test data of clean and artificially iced aircraft, and a modified stepwise regression algorithm was used to obtain pitching and yawing derivatives. Whalen [6], based on the flight test data of Twin Otter, studied the changes of

Aircraft Dynamic Model
Parameter identification is to solve the related parameters of a system when the input and output of the system are known. To describe the system, a mathematical model is built first. According to the characteristics of an iced aircraft, a nonlinear longitudinal flight dynamics model is constructed.
where the state vector → x (t) includes forward velocity u, angular velocity w, q, the Euler angle θ and spatial position coordinates x e , h. The control vector → u (t) includes the elevator deflection δ e and engine thrust δ T . Figure 1 shows the coordinate systems of longitudinal motion of airplane, where ox b and oz b are the body axes, ox E and oz E are the earth axes.
x e h = L Eb u w (2) where L Eb is the transfer matrix from the body coordinate system to earth's coordinate system.  e Eb where Eb L is the transfer matrix from the body coordinate system to earth's coordinate system. Generally, aircraft is regarded as a rigid body when analyzing the whole motion of aircraft. The nonlinear ordinary differential motion equations are as follows: The aerodynamic model is as following: Generally, aircraft is regarded as a rigid body when analyzing the whole motion of aircraft. The nonlinear ordinary differential motion equations are as follows: The aerodynamic model is as following: where, C X0 , C Xα2 , C Z0 , C Zα , C Zq , C Zδe , C m0 , C mα , C mq , C mδe are unknown aerodynamic derivatives, that is, the target of parameter identification.

Parameter Identification Method
After the aircraft dynamic model is established, the unknown parameters of the model, namely aerodynamic derivative, need to be identified. In this paper, the maximum likelihood criterion is used to identify the parameters.
The dynamic system model can be written as follows: The theory of parameter identification proves that the maximum likelihood estimation of parameters is asymptotically unbiased, consistent and effective. When the experimental data of the system are enough, the mathematical expectation of the parameters obtained by the maximum likelihood identification is the true value and converges to the true value with probability 1. Its variance gradually reaches the Cramer-Rao lower bound, which is the unbiased estimation with the minimum variance. Therefore, the maximum likelihood estimation of the parameters is to seek the parameters to minimize the following likelihood indicator functions.Indicator function of the maximum likelihood identification is: where B is the innovation variance matrix. The optimal estimationB of B can be obtained by finding the extremum of J to B.B Then, the asymptotically unbiased, uniform and efficient estimation of → λ can be obtained by minimizing J to → λ. The iterative convergence process of → λ is based on the modified Newton-Raphon algorithm [19].
The procedure of the identification program is shown in Figure 2.

The Prediction Method
The maximum likelihood identification method is introduced to identify the mathematical model under aircraft icing in Section 2.2. However, this method strongly relies on the flight test, whenever the ice shapes have a large change the identification results maybe heavily discounted. Therefore, based on the values of clean aircraft's longitudinal aerodynamic derivatives and the ice shapes of wing and horizontal tail, a prediction method is established to calculate the longitudinal aerodynamic derivatives of C Zα , C Zq , C Zδe , C mα , C mq and C mδe . The correctness of the prediction method will be verified by comparing the calculated results with the identified results. Figure 2. Procedure of the identification program.

The Prediction Method
The maximum likelihood identification method is introduced to identify the mathematical model under aircraft icing in Section 2.2. However, this method strongly relies on the flight test, whenever the ice shapes have a large change the identification results maybe heavily discounted. Therefore, based on the values of clean aircraft's longitudinal aerodynamic derivatives and the ice shapes of wing and horizontal tail, a prediction method is established to calculate the longitudinal aerodynamic derivatives of Z C  , Zq C , Z e C  , m C  , mq C and m e C  . The correctness of the prediction method will be verified by comparing the calculated results with the identified results. The aerodynamic force in z direction on clean aircraft could be expressed as following: is the interference aerodynamic force caused by other part of aircraft components.
The aerodynamic force in Z direction on aircraft with iced wing and horizontal tail could be expressed as following: Assume that the changes of aerodynamic force on body and the mutual influence of aircraft components caused by icing could be ignored, that is: , , The aerodynamic force in z direction on clean aircraft could be expressed as following: where Z (inter) is the interference aerodynamic force caused by other part of aircraft components. The aerodynamic force in Z direction on aircraft with iced wing and horizontal tail could be expressed as following: Assume that the changes of aerodynamic force on body and the mutual influence of aircraft components caused by icing could be ignored, that is: Then, Equation (11) can be obtained from Equations (8)- (10).
namely, (12) where k q(ice) and k q represent the ratios between the dynamic pressure on the horizontal tail and the wing, and the estimation method from [20] is used to calculate them. Obtain the partial derivatives of Equation (12) with respect to α, δe and q, separately, and ignore the effects of them on k q(ice) and k q . Then Equation (13) could be achieved.
Assume that the changes of C Zα , C Zδe and C Zq caused by the change of q ∞ could be ignored in a certain range. Suppose that In addition, because k q is the ratios of dynamic pressure q ∞ and it is assumed that the dynamic pressure q ∞ is constant before and after icing. The effects of icing on k q could be ignored, namely, Substituting Equations (14) and (15) into (13), Equation (16) could be obtained.
In a similar manner to the upper derivation, Equation (17) could be obtained.
Ignore the change of aerodynamic force location on wing and horizontal tail caused by icing. Equation (18) could be obtained from Equation (17).
In Equations (16) and (18), to obtain the difference between iced and clean aircraft aerodynamic derivatives, the terms in the right hand of the equation should be solved first. C Zα , C Zδe could be easily calculated by changing the angle of attack or the elevator location using CFD software.
As for the dynamic derivative of C Zq , a calculation method called the 'strip' technique developed by [21] is used. The formulas are shown as Equation (19).

Identification Method Accuracy Verification
In order to verify the accuracy of the identification method, the following steps are taken.
(1) Choose a set of aerodynamic parameters → λ f as the reference value under a flight condition within a reasonable range.
(2) Take → λ f as the parameters of an aircraft dynamic model to obtain the response of this reference dynamic model under certain elevator inputs. Then, regard the response data as the flight test reference response for identification using in the next step.
Start this identification process from the initial value → λ 0 using the response data obtained in step 2. The input value → λ r is calculated iteratively until the input value response is basically consistent with the reference response. Then, the result → λ r can be achieved.
(4) By comparing the difference between the identification result and the reference value, the correctness and accuracy of the identification system are verified.
Aircraft DHC-6 Twin Otter [18], a high-wing, twin-engine, commuter-class aircraft, is chosen as the test plane. The initial flight conditions are shown in Table 1. The result of aerodynamic parameters is shown in Table 2, and the output responses fitting curve is shown in Figure 3, respectively. As shown in Table 2, the results compared extremely well and most of the errors between identification results and reference values are less than 6%. The errors of Zq C and Z e C  are larger than others, 12.50% and 11.17%, respectively. However, the change tendency is reasonable.  As shown in Table 2, the results compared extremely well and most of the errors between identification results and reference values are less than 6%. The errors of C Zq and C Zδe are larger than others, 12.50% and 11.17%, respectively. However, the change tendency is reasonable.
It could be observed from Figure 3 that the initial responses of angle of attack and pitch rate have some discrepancies as compared with the reference response curve. However, when the system parameters are alternated with the identification result → λ r , the resulting response curves are almost coincident with the reference curves.
According to the verification results, the error of the identification method is within the acceptable range. This method is relatively effective and accurate, and can be used to identify the aircraft longitudinal aerodynamic parameters.

Identification of Clean Aircraft
Firstly, the identification of clean aircraft in doublet pair flight is carried out. The results are used as the basis for analyzing the effect of aircraft icing. The aerodynamic parameters of aircraft DHC-6 are achieved according to the flight test data provided by [6] (Doublet Pair Flight). The initial flight conditions of clean aircraft are shown in Table 3 and the output responses fitting curves are shown in Figure 4. Figure 4 shows the responses of the system model obtained by the identification accord quite well with the flight test data. This proves that the aircraft dynamic model established in this paper can simulate DHC-6 aircraft relatively accurately. Furthermore, it proves that this identification method is useful for identifying aerodynamic parameters.
The initial value and identification results are compared in Table 4. The results reported in [10] are taken as the references.  It could be observed from Figure 3 that the initial responses of angle of attack and pitch rate have some discrepancies as compared with the reference response curve. However, when the system parameters are alternated with the identification result r   , the resulting response curves are almost coincident with the reference curves. According to the verification results, the error of the identification method is within the acceptable range. This method is relatively effective and accurate, and can be used to identify the aircraft longitudinal aerodynamic parameters.

Identification of Clean Aircraft
Firstly, the identification of clean aircraft in doublet pair flight is carried out. The results are used as the basis for analyzing the effect of aircraft icing. The aerodynamic parameters of aircraft DHC-6 are achieved according to the flight test data provided by [6] (Doublet Pair Flight). The initial flight conditions of clean aircraft are shown in Table 3 and the output responses fitting curves are shown in Figure 4.     Figure 4 shows the responses of the system model obtained by the identification accord quite well with the flight test data. This proves that the aircraft dynamic model established in this paper can simulate DHC-6 aircraft relatively accurately. Furthermore, it proves that this identification method is useful for identifying aerodynamic parameters.
The initial value and identification results are compared in Table 4. The results reported in [10] are taken as the references.   As shown in Table 4, although the identification results have some discrepancies relative to the initial values, most of the aerodynamic derivatives accord with the reference results quite well.
The causes of the discrepancies may be as follows: (1) In different flight tests, the weight and thrust coefficient of DHC-6 aircraft are different.
(2) In the process of identification, the aircraft dynamic model established in this paper is different from the reference.

Identification of Horizontal Tail Iced Aircraft
The horizontal tail is the important part of aircraft for the ability to affect aircraft longitudinal control and stability. To study the effect of icing on horizontal tail aerodynamic performance, data from [18] are used, and the ice shapes of "Inter-cycle" and "Stability and Control (S&C)" on the horizontal tail are chosen, as shown in Figure 5. The identification results of clean aircraft in Section 4.1 are used as the initial value. As shown in Table 4, although the identification results have some discrepancies relative to the initial values, most of the aerodynamic derivatives accord with the reference results quite well.
The causes of the discrepancies may be as follows: (1) In different flight tests, the weight and thrust coefficient of DHC-6 aircraft are different.
(2) In the process of identification, the aircraft dynamic model established in this paper is different from the reference.

Identification of Horizontal Tail Iced Aircraft
The horizontal tail is the important part of aircraft for the ability to affect aircraft longitudinal control and stability. To study the effect of icing on horizontal tail aerodynamic performance, data from [18] are used, and the ice shapes of "Inter-cycle" and "Stability and Control (S&C)" on the horizontal tail are chosen, as shown in Figure 5. The identification results of clean aircraft in Section 4.1 are used as the initial value. The S&C ice shape is used extensively in NASA's previous stability and control tests and it is derived from in-flight photos and ADS-4. The aeroperformance characteristics for this ice shape are the worst and nearly identical to the "lewis-ice (LEWICE)" shape. LEWICE shape conditions: V = 120 knots, alpha = −2.9°, LWC = 0.5 g/m 3 , MVD = 20 um, T0 =−4°C, time = 45 min [18]. Table 5 shows the initial flight conditions of horizontal tail iced aircraft. The fitting curves of the output responses of the aircraft with two different ice shapes are shown in Figure 6, and the result of aerodynamic parameters are shown in Table 6, respectively. The S&C ice shape is used extensively in NASA's previous stability and control tests and it is derived from in-flight photos and ADS-4. The aeroperformance characteristics for this ice shape are the worst and nearly identical to the "lewis-ice (LEWICE)" shape. LEWICE shape conditions: V = 120 knots, alpha = −2.9 • , LWC = 0.5 g/m 3 , MVD = 20 um, T 0 =−4 • C, time = 45 min [18]. Table 5 shows the initial flight conditions of horizontal tail iced aircraft. The fitting curves of the output responses of the aircraft with two different ice shapes are shown in Figure 6, and the result of aerodynamic parameters are shown in Table 6, respectively.
According to the identification results of the two kind of horizontal tail iced aircraft as shown in Table 6, the following analysis could be made.
(1) Compared with clean aircraft, the absolute value of C X0 increases by 27.03% and 33.24%, respectively, under the condition of horizontal tail with two ice shapes. This means that aircraft icing will greatly increase aircraft drag. The absolute values of C Z0 and C Zα are slightly reduced, which indicates that horizontal tail icing will lead to a small reduction of lift. (2) The absolute value of C mq changes largely under icing conditions with Inter-cycle ice shape and S&C ice shape, by −14.62% and −20.27%, respectively, which means the horizontal tail icing would cause the decrease of pitching damping. (3) Relative to the clean aircraft, the absolute values of C Zδe and C mδe for S&C ice shape condition change significantly by −46.92% and −21.29% respectively, which indicates that horizontal tail icing contributes significantly to the reduction of elevator effectiveness.
(4) Comparing the identification results under the two different icing conditions, the aerodynamic parameters with S&C ice shape changes larger than those with Inter-cycle ice shape, which means horizontal tail iced aircraft with the ice shape of S&C is more affected.   (c) Elevator deflection 2 [18] (d) S&C ice shape

Identification of All Configuration Iced Aircraft
Identification of all configuration iced aircraft is studied in this section. With the aerodynamic parameters of clean aircraft in Section 4.1 as the initial values and the flight test operation (Elevator Doublet) in [17] as the input values, the aerodynamic parameters of "all iced" aircraft are identified. These flight tests were conducted with the whole iced configuration of Twin Otter's aircraft, including the wing, horizontal tail, and vertical tail. The ice shapes of those plans could be found from [17] as shown in Figure 7.
would cause the decrease of pitching damping.
(3) Relative to the clean aircraft, the absolute values of Z e C  and m e C  for S&C ice shape condition change significantly by −46.92% and −21.29% respectively, which indicates that horizontal tail icing contributes significantly to the reduction of elevator effectiveness. (4) Comparing the identification results under the two different icing conditions, the aerodynamic parameters with S&C ice shape changes larger than those with Inter-cycle ice shape, which means horizontal tail iced aircraft with the ice shape of S&C is more affected.

Identification of All Configuration Iced Aircraft
Identification of all configuration iced aircraft is studied in this section. With the aerodynamic parameters of clean aircraft in Section 4.1 as the initial values and the flight test operation (Elevator Doublet) in [17] as the input values, the aerodynamic parameters of "all iced" aircraft are identified. These flight tests were conducted with the whole iced configuration of Twin Otter's aircraft, including the wing, horizontal tail, and vertical tail. The ice shapes of those plans could be found from [17] as shown in Figure 7. (a) Ice shape of wing (b) Ice shape of horizontal tail Figure 7. Ice shape of "all iced" aircraft [17].
The initial flight conditions of "all iced" aircraft are shown in Table 7. This initial flight conditions are quite close to the conditions of clean aircraft, which will help to get a better identification result. The output responses fitting curves are shown in Figure 8, and the identification results are shown in Table 8, respectively. The initial flight conditions of "all iced" aircraft are shown in Table 7. This initial flight conditions are quite close to the conditions of clean aircraft, which will help to get a better identification result. The output responses fitting curves are shown in Figure 8, and the identification results are shown in Table 8, respectively.  Value 1200 m 55 m/s 0.14 (a) Elevator deflection [17] (b) Response of angle of attack   According to the identification results of artificially iced aircraft as shown in Table 8, the following analysis could be made.

Changes of Horizontal
(1) Under an "all iced" condition, the absolute value of C X0 increased by 67.57%. The absolute values of C Z0 and C Zα decreased by −8.44% and −6.13%, respectively. Compared with the horizontal tail icing condition, these parameters change more. The change range is more than twice. The results show that the main wing icing plays a leading role in the change of aircraft lift and drag characteristics under icing conditions. (2) The changes of C Zq and C mq do not reflect great difference compared with horizontal tail iced aircraft, which means when the airplane has a tail, the wing contribute to C Zq and C mq is often negligible in comparison with that of tail [22]. (3) C mα changes largely under "all iced" condition by −30.29%, while those of horizontal tail iced condition with S&C and Inter-cycle ice shape are −4.15% and −11.34%, respectively. These results indicate that wing icing contributes more to the change of C mα . (4) The changes of C Zδe and C mδe under the "all iced" condition are slightly more than those under a horizontal tail iced condition, which indicate that horizontal tail icing mainly contributes to the reduction of elevator effectiveness.
Although the results of this identification method are well fitted with the experimental results and have high accuracy, there are some limitations. This identification method is highly dependent on flight test. When the aircraft configuration is different from that of the identification model, or the ice shape is greatly changed, the accuracy of the identification result will be greatly reduced. Therefore, it is necessary to propose an engineering prediction method to calculate the aerodynamic derivatives of the aircraft without sufficient flight test data.

Prediction Method Application and Result Analysis
The aerodynamic parameters of the wing and horizontal tail: C ZαW , C ZαH , C ZδeH are calculated under both clean and iced conditions with the ice shapes of the wings and horizontal tail chosen from [17] by the individual component CFD method. As an illustration, Figure 9 shows (2) The changes of Zq C and mq C do not reflect great difference compared with horizontal tail iced aircraft, which means when the airplane has a tail, the wing contribute to Zq C and mq C is often negligible in comparison with that of tail [22]. (3) m C  changes largely under "all iced" condition by −30.29%, while those of horizontal tail iced condition with S&C and Inter-cycle ice shape are −4.15% and −11.34%, respectively. These results indicate that wing icing contributes more to the change of m C  .
(4) The changes of Z e C  and m e C  under the "all iced" condition are slightly more than those under a horizontal tail iced condition, which indicate that horizontal tail icing mainly contributes to the reduction of elevator effectiveness.
Although the results of this identification method are well fitted with the experimental results and have high accuracy, there are some limitations. This identification method is highly dependent on flight test. When the aircraft configuration is different from that of the identification model, or the ice shape is greatly changed, the accuracy of the identification result will be greatly reduced. Therefore, it is necessary to propose an engineering prediction method to calculate the aerodynamic derivatives of the aircraft without sufficient flight test data.

Prediction Method Application and Result Analysis
The aerodynamic parameters of the wing and horizontal tail:  According to the calculation results, the curves of ZH C with angle of attack are drawn, as shown in Figure 10. The slope obtained by linear fitting the two curves can be used as the estimated value of ZH C . The other parameters are calculated in the same way. According to the calculation results, the curves of C ZH with angle of attack are drawn, as shown in Figure 10 The results are shown in Table 9. Additionally, the values of zqW C and zqH C can be calculated from Equation (19). The value of differences between iced and clean aircraft aerodynamic parameters can be obtained by replacing the corresponding terms in Equations (16) and (18) with the above-calculated results. Then the aerodynamic parameters of iced aircraft are obtained by adding the differences to the aerodynamic parameters of clean aircraft identified in Section 4.1. The results are compared with the identification results of iced aircraft in Section 4.2 as shown in Table 10.  The results are shown in Table 9. Additionally, the values of C zqW and C zqH can be calculated from Equation (19). The value of differences between iced and clean aircraft aerodynamic parameters can be obtained by replacing the corresponding terms in Equations (16) and (18) with the above-calculated results. Then the aerodynamic parameters of iced aircraft are obtained by adding the differences to the aerodynamic parameters of clean aircraft identified in Section 4.1. The results are compared with the identification results of iced aircraft in Section 4.2 as shown in Table 10. As shown in Table 10, the change tendency of the prediction results is coincident with the identification results. For several aerodynamic derivatives of iced aircraft, the error between calculation results and identification results is basically within 10%. Only C Zδe reaches 15.84%. These errors may be caused by some assumptions or simplifications used in the prediction method. For the dynamic derivative C mq , which represents the pitch damping, the error 3.55% is not very large. Thus, the feasibility of the prediction method in estimating the longitudinal dynamic derivatives of iced aircraft is verified. Based on these validations, the prediction method with the coupled 'strip' technique and CFD calculation shows good agreement with the experiment, and it could be used to calculate the effects of icing on the aircraft longitudinal characteristics. With this prediction method, only the aircraft basic configuration, the clean aircraft longitudinal aerodynamic derivatives and the ice shape are needed to calculate the longitudinal aerodynamic derivatives of the iced aircraft. Therefore, this method has a wide range of application.

Conclusions
In this paper, a nonlinear longitudinal flight dynamics model of rigid aircraft is established. For the aircraft with flight test, based on the maximum likelihood criterion, the identification system of aircraft aerodynamic parameters is established. For the aircraft without a flight test, based on the individual component CFD calculation and narrow strip theory, an engineering prediction method of longitudinal aerodynamic derivatives is established. According to NASA's icing flight test data, longitudinal aerodynamic parameters of clean and icing aircraft in different conditions were identified. The aerodynamic parameters of the same icing aircraft are also calculated by the engineering prediction method. The difference between the results of the two methods and the reference value is compared, and the correctness of the two methods is verified. Through the above research and analysis, we can draw the following conclusions: (1) The aerodynamic parameters identification system based on the maximum likelihood criterion can predict the longitudinal aerodynamic parameters of iced aircraft relatively accurately. Through the analysis and verification of several examples, the method is practical. (2) The engineering prediction method established in this paper shows good agreement with the experiment. On the basis of known aircraft configuration and ice shape, this method can be used to calculate the effect of icing on the aircraft longitudinal aerodynamic parameters.