#### Abstract

Current fight boundary of the envelope protection in icing conditions is usually defined by the critical values of state parameters; however, such method does not take the interrelationship of each parameter and the effect of the external disturbance into consideration. This paper proposes constructing the stability boundary of the aircraft in icing conditions through analyzing the region of attraction (ROA) around the equilibrium point. Nonlinear icing effect model is proposed according to existing wind tunnel test results. On this basis, the iced polynomial short period model can be deduced further to obtain the stability boundary under icing conditions using ROA analysis. Simulation results for a series of icing severity demonstrate that, regardless of the icing severity, the boundary of the calculated ROA can be treated as an estimation of the stability boundary around an equilibrium point. The proposed methodology is believed to be a promising way for ROA analysis and stability boundary construction of the aircraft in icing conditions, and it will provide theoretical support for multiple boundary protection of icing tolerant flight.

#### 1. Introduction

Aircraft icing has always been a familiar but unsolved factor threatening flight safety. The formation of ice on aircraft in flight alters the original smooth aerodynamic shape of the aircraft, which has a negative effect on aircraft performance, stability, and control. The shrinking flight envelope caused by aircraft icing can easily lead to flight accident. Obviously, ice avoidance is the most reliable way to increase flight safety; however this conservative way is unnecessary. To improve the operation ability of the aircraft and to maintain revenues and schedules, ice tolerance is the preferred choice for icing conditions.

The goal of the ice tolerance flight is to make sure that the aircraft operate in the safe flight envelope. For this purpose, great efforts have been made by many researchers around the world. Bragg et al. proposed a Smart Icing System conception [1, 2]; the system could sense and characterize the presence of ice, notify the pilot, and ensure the safety of the aircraft. Sharma et al. researched the envelope protection theory when operating in autopilot mode as one part of SIS research [3]. Gingras et al. [4, 5] have developed the Icing Contamination Envelope Protection (ICEPro) system which is mainly designed to provide flight and control limit information to the pilot through cues and messaging on his display to improve flight safety. All above achievements provide systematic approaches for ice tolerance flight. In addition, researches related to one aspect of ice tolerance have also been carried out to perfect the theory and methodology.

Melody et al. [6] compared three different parameter identification algorithms in the context of icing detection and pointed out that method can provide a timely and accurate icing indication. Caliskan et al. [7] adopted Kalman filter and neural network technique to study problems on aircraft icing identification, detection, and reconfigurable control. Dong and Ai proposed using time-varying case of the algorithm and probabilistic neural network to provide inflight parameter identification and icing location detection [8].

The definition of the flight boundaries under icing conditions is one of the key aspects of ice tolerance flight. However, current flight boundaries are mainly defined by critical values of several state parameters. For example, the decreased stall angle of attack due to icing is commonly used to restrict the incidence angle within its safe envelope. Such method has two main limitations: on the one hand, no matter single parameter limit or several parameters limit, the interrelation of each key parameter is not taken into consideration, while loss of control (LOC) may happen caused by the coupling of state parameters although each of them is within its own limitation. One the other hand, the method is feasible when neglecting the outside disturbance, whereas the aircraft can always be perturbed due to an upset condition or wind gust which may lead to system instability even if the state parameter does not exceed its extreme value. To keep the passengers comfort and safety flight under icing conditions, it is of great importance to make sure that the aircraft operates in the stable flight domain.

The paper presented here aims at developing a methodology to conduct a preliminary analysis of the stability boundary of the aircraft under icing conditions. The remainder of the paper is organized as follows. Section 2 proposes a nonlinear icing effect model for a wide range of angle of attack to incorporate the icing effect into flight dynamics. The process of derivation of the longitudinal axis polynomial model of the aircraft is presented in Section 3. Section 4 describes the underlying theory of region of attraction (ROA) analysis for nonlinear polynomial system based on sum-of-square (SOS) theory. Comparison between ROA analysis results and the phase diagram of different icing severity presented in Section 5 demonstrate that the boundary of the calculated ROA can be regarded as the approximate flight stability boundary of the aircraft under icing conditions. Concluding remarks and recommendations for the future work are presented at the end of this paper.

#### 2. Nonlinear Icing Effect Model

Modeling of abnormal aerodynamics caused by ice accretion has been studied for many years [9–12]. Most of these researches only focus on linear interval of incidence angle, that is, within stall angle. To give an insight into flight dynamic characteristics near the flight boundary especially for poststall region, nonlinear icing effect model for a wide range of angle of attack should be established. In this section, a nonlinear icing effect model is proposed to include the adverse effect of ice accretion on aerodynamic parameters.

##### 2.1. Linear Icing Effect Mode

Current estimation of iced aerodynamics model normally based on linear stability and control derivatives, namely, aerodynamic force or moment, is linear function of flight state parameters and control value. For longitude aerodynamic parameters, the aerodynamic force coefficients acting along the body axes, , and the pitch moment coefficients are calculated by

Lift coefficient and drag coefficient can be available via a rotation

For linear icing effect model, the effect of ice accretion is reflected by the change of the aerodynamic derivatives. The icing effect model proposed by Bragg et al., which combines aircraft specific property and icing conditions together, is physically representative and has been widely used in icing research [12]. According to the theory, the icing effect theory is based on the following equation:where is any arbitrary aerodynamic coefficient that is affected by ice accretion. represents icing severity which is only the function of the atmospheric conditions. is a constant value which denotes the modification of the aerodynamic coefficients; it is dependent on the coefficient being modified and on the properties of the specific aircraft.

##### 2.2. Nonlinear Icing Effect Model

###### 2.2.1. Establishment of the Model

The adopted equation form of the longitude nonlinear aerodynamics model is [2, 13]where , , and are unitless aerodynamic coefficients computed from look-up tables.

Compared with the linear aerodynamic model, both of the two models have similar formation. For example, as for , the term in linear aerodynamic model (see (1)) is corresponding to the term in nonlinear aerodynamic model (see (1)). And , terms are corresponding to terms, respectively. According to linear icing effect model shown as (3), each term of in the nonlinear aerodynamics model can be calculated bywhere subscription iced and clean denotes the flight condition with and without ice, respectively.

###### 2.2.2. Modification of Nonlinear Icing Effect Model

According to the subscale, complete airplane wind tunnel test for a wide range of angle of attack, the differences between clean and iced configuration are gradually narrowing and finally disappear with increasing angle of attack in poststall area for longitudinal aerodynamic coefficients [14–17]. Typical alteration characteristics especially in high angle of attack are summarized as follows:(i)A delayed nose-down break for pitch moment coefficient in the stall region.(ii)A slight reduction in lift curve slope and a sizable reduction in maximum lift.(iii)A “flattening” of the lift curve slope in the stall region.(iv)Lift coefficient decrease or increase again at almost linear lift curve slope in the poststall region.

However, the direct calculated results through (11) shown by blue line in Figure 1 cannot reflect those phenomena, and the disagreement is due to the invariant value of . In fact, is no longer a constant value when it is in high angle of attack and it is also the function of state parameters like [2]. Based on above analysis, modification on iced coefficients can be divided into following intervals:(1)For linear change interval, is the same as linear icing effect model.(2)In the stall region, value of is depend on the corresponding aerodynamic coefficient, such as curvature of lift coefficient curve in stall region decrease as the severity of ice increase.(3)In poststall region, the variation of for iced aerodynamic coefficient is determined by the approximate linear relationship of aerodynamic coefficient and angle of attack. When the angle of attack reaches the certain value, turn out to be zero, which means the differences between clean and iced configuration disappear through point.

**(a)**

**(b)**

**(c)**As to drag coefficient term , which is calculated through (2) and (5) and without modification, the computed results shown in Figure 1(b) indicate that the drag coefficient decreases when ice accretion occurs in certain interval. This is quite unreasonable for ice accretion which generally leads to drag increase. To obtain the rational variation trend, the calculation of should be counted as separate problem. The computing method is also based on (3) but the value of depends on the angle of attack.

The estimation of icing severity parameter is not within the research scope of this paper, and the detailed process can refer to [12].

Most of the current icing effect model can only get the decreased lift coefficient and cannot reflect the change of stall angle. This is due to the essence of these models; they only count the percentage of alteration for each derivative. Although the value of the maximum lift coefficient comes out to be reduced, the stall angle is still the same. The commonly used solution is calculating the stall angle separately when, in flight simulation, this method is feasible when . However, when calculating the flight dynamic characteristics in stall and poststall area, the method is not proper. To solve this problem, the reduction factor is defined to measure the alteration of the stall angle. By the definition of , the original incidence angle interval can be mapped into a new interval as follows:where and it has a negative correlation with icing severity. is the initial point in the look-up table and is the set point to determine the interval to be shrink. In this way, the reduced stall angle can be incorporated into icing effect model. Adversely, the delayed nose-down break for pitch moment coefficient can be obtained by extending the specific angle interval. where and it has a positive correlation with icing severity.

Based on the proposed nonlinear icing effect model, comparison between clean and iced aerodynamic coefficients with and without modification for a typical transportation aircraft is shown in Figure 1. For brevity, only static parameters are shown.

In Figure 1, the black line shows the clean aerodynamic coefficient and the red line shows the modified iced aerodynamic coefficient. Compared with the wind tunnel data of airfoil and subscale, complete airplane models [14–17], the established nonlinear icing effect model has similar variation tendency with the experiment results for a wide range of angle of attack.

#### 3. Polynomial Aircraft Model under Icing Conditions

##### 3.1. Conventional Longitudinal Dynamic Model

The conventional longitudinal nonlinear dynamic model of an aircraft can be described as follows [18]:where* m* is mass,* U* is the airspeed, is the angle of attack, is the pitch angle, and is the pitch rate. and are projection of the thrust along the body axes and , respectively, and denotes the pitch moment produced by thrust which is calculated by function of and . The lift force* L*, drag force* D*, and pitch moment* M* in the model are given bywhere is the dynamic pressure and is the normalized pitch rate (unitless). is unitless aerodynamic coefficient and denotes the mean aerodynamic chord.

###### 3.1.1. Polynomial Longitudinal Model

As shown in the above conventional nonlinear dynamic model in (9) and (10), there are several nonpolynomial terms like trigonometric functions term, derivative airspeed item, and aerodynamic coefficient term that could not be directly used to ROA analysis based on sum-of-square theory. Therefore, it is necessary to transfer these nonpolynomial terms into polynomial terms before carrying on ROA analysis. It is obvious that the accuracy of the transferred polynomial increases with the polynomial order. However, the calculation of the ROA will take more time either. Thus, the proper order of the transferred polynomial is very important.

The trigonometric functions can be converted into polynomial through third-order Taylor series expansion:

Error is in a reasonable range for and the maximum relative error in sine and cosine is 0.54% and 3.08%, respectively, which demonstrate that the transferred polynomial is accurate enough for substituting the trigonometric function term.

For derivative airspeed item, 2nd-order polynomial fit is used to obtain the approximate polynomial expression (12) on the specified interval.where* a*,* b*, and* c* are constant coefficients to be determined.

It is assume that the engine thrust is only relevant to the throttle position ; thus the polynomial model of engine thrust can also obtained by polynomial fit. Order of the model depends on the specific function relationship between thrust value and throttle position, which is generally no more than three.

As shown in (10) aerodynamic coefficient is commonly related to flight state parameters (like* α* and

*q*) and control variables (like ). Therefore, multivariate polynomial fit is needed. For iced aerodynamic coefficient curves, as shown in Figure 1, are more twist and turn, higher degree of polynomials is required to fit the iced coefficients better. For example, polynomial model of , , and can be given bywhere are constant coefficient to be determined. To ensure trim characteristics of the aircraft, weighted least square fit method [18] is adopted to calculate the coefficient.

After substituting all the nonpolynomial terms in conventional longitudinal model (see (9)), polynomial model can be available as follows:where and .

###### 3.1.2. Modeling and Verification of Iced Polynomial Model

For illustrating, Generic Transport Model (GTM) is used to verify the effectiveness of the presented icing effect model and to research the variation of the stability boundary in different icing conditions.

Feasibility of the polynomial model has been verified in the previous research [13, 18]. In this section, the trim conditions and dynamic characteristics of the polynomial model that incorporated with proposed icing effect model are computed to assess the quality of polynomial model and the icing effect model. The trim conditions are assumed to keep level flight at* H *= 500 m and* V *= 50 m/s. Figure 2 shows the trim angle of attack and trim inputs versus different icing severity for both the conventional nonlinear model and the polynomial model. Consistency of the trim behavior of these two models demonstrates the effectiveness of the polynomial model in aircraft icing research to certain extent.

Figure 2 also shows that the trim angle of attack and trim inputs increase as ice accretion became more severe. The phenomenon is due to the increased drag, reduced lift coefficient, and control effectiveness caused by aircraft icing. This is consistent with previous research, which demonstrates the effectiveness of the proposed icing effect model in some extent.

After the aircraft is being perturbed due to a wind gust or other upset conditions, parameters of the short period mode change rapidly and once they divergent, they will make it hard to ensure flight safety. On the other hand, parameters of the phugoid mode change slowly enough for the pilot to modify the oscillation of the aircraft. Generally, when conducting stability analysis, only short period mode parameters are considered for simplicity.

For polynomial model, the short period model can be extracted from the 4-state model by holding* V* and* θ* and at their trim values , , and :where are the second and third formulations of the conventional 4-state model, (14), respectively. For the clean GTM model, the corresponding polynomial short period model comes out to be

In this section, simulations in different initial conditions of polynomial short period model and the conventional short period model are performed to further validate the effectiveness of using polynomial short period model in ROA analysis. Initial condition of the simulations of the two model are set to be the same with the elevator inputs held fixed at its trim value. For comparison, trajectories of parameters* α* and

*q*are shown on - phase plane. As the valid range of angle of attack is larger than −5°, only situations are demonstrated here as shown in Figure 3.

The calculated result shows a good agreement between these two models when flight state parameters stay in their valid range. Thus, it is feasible to use the polynomial short period model to research the stability characteristic of the GTM model.

According to the above validation procedure, as long as the flight parameters* α* and

*q*are in the valid scope of the polynomial model, the polynomial short period model can adequately capture the short period stability characteristic of the origin model. As the number of the state parameters and the degree of the polynomial model will greatly affect the computational complexity, the polynomial short period model is preferred when conducting ROA analysis.

#### 4. ROA Estimation Based on SOS Theory for Nonlinear Polynomial System

Region of attraction is an important tool in nonlinear dynamics analysis and it has been used for flight control analysis, envelope protection, and so forth [18–21]. ROA of a locally asymptotically stable equilibrium point is an invariant set such that all trajectories starting inside this set converge to the equilibrium point [22]. For an aircraft, the size of the ROA means that the extent of flight states deviates away from the equilibrium point, which can be driven by wind gust or other upset conditions, and instability does not occur. As long as the flight state is in the area of ROA, the airplane will converge back to its initial trimmed flight state. The smaller the ROA is, the more easier the system is affected by outside disturbance. Therefore, taking the boundary of ROA into envelope protection as limitation of state parameters will bring a more reliable safety flight area.

Consider the autonomous nonlinear system [22]where is the state vector and is a vector polynomial function of with . Assume that the origin is a locally asymptotically stable equilibrium point; the ROA of the origin is defined aswhere is the solution of (1) that starts at initial state .

Generally, the exact ROA for nonlinear dynamical systems is very hard to compute and current research can only determine the invariant subsets of the ROA [22, 23]. In this paper, Lyapunov function based method is adopted. This kind of method computes a Lyapunov function (LF) as a local stability certificate and sublevel sets of this LF provide invariant subsets of the ROA [24]. Numerical algorithm in this section is based on Lemma 1 that will be described in the following.

Lemma 1. *If there exists and a polynomial such that [18]*(1)* and ;*(2)* is bounded;*(3)*,**then for all , the solution of (17) exists and . As such, is a subset of the ROA for (17). A function V satisfying the conditions in Lemma 1 is a Lyapunov function and provides an estimate of the region of attraction. To enlarge , a variable sized region is defined as follows:and it maximizes under the constraint . is a positive value and is shape function which is normally defined as . The choice of reflects the relative importance of states that may have great influence on the computed ROA [22]. According to Positivstellensatz theorem [25], the above optimization problem can be transferred and further simplified as following Sum of Square Programming (SOSP) problem:where is the set of all polynomials in n variables and is the set of SOS polynomials in n variables. are positive definite polynomial of the formfor and are small positive constants on the order of . In this paper, - iteration algorithm [26] is adopted to solve the SOSP. Finally, the stability boundary of the polynomial nonlinear system can be depicted as *

*Apparently, as long as the flight state inside the area which is enclosed by the stability boundary, the aircraft will converge to the trim condition.*

#### 5. Stability Boundary in Different Icing Conditions

##### 5.1. Representation of the ROA in Different Icing Conditions

After obtaining the polynomial short period model of the aircraft, ROA and stability boundary can be worked out based on ROA estimation theory in Section 2. To research the variation of the stability boundary when the aircraft encounters ice accretion, ROA analysis is conducted for different icing conditions in this section. For all cases, the trim condition is set to be altitude = 500 m and velocity = 50 m/s.

For brevity, two modifications of the polynomial model are adopted. The first one is shifting the nonzero trim state to the origin point to make the polynomial model suitable for ROA analysis algorithm. The second one is removing all polynomial terms with degree greater than five and coefficients less than 10^{−8} to reduce the computation time [18].

For clean configuration, the modified polynomial short period model can be worked out as follows:

Based on SOS theory, ROA around the origin point is

Thus the estimated stability boundary is

Similarly, the modified polynomial short period model and the corresponding ROA in different icing conditions , 0.2, and 0.3 are also computed; detailed results are provided in Appendix. To further illustrate the variation of the stability boundary in different icing conditions, all of the above estimated stability boundaries are shown in Figure 4. As discussed in Section 4 the calculated ROA is only valid in the range of due to the limitation of the data available. Thus, the effective region of the ROA is enclosed by vertical line (shown in green dashed line) and right part of as shown in Figure 4. Apparently, as the severity of aircraft icing increases, the area, enclosed by the estimated stability boundary, becomes smaller; that is, the aircraft is more apt to loss stability induced by outside disturbance.

Take the effective area of ROA as an indicator to measure the variation of the boundary. For clean configuration, the effective area of the ROA worked out as 7.2692 × 10^{3}, for iced configuration of , 0.2, and 0.3, and the effective area of the ROA came out to be 5.5576 × 10^{3}, 3.6134 × 10^{3}, and 2.1872 × 10^{3}, respectively (neglect the unit). Further, as the ice accretion deteriorates, the size of the estimated stability region reduces by 23.54%, 50.29%, and 69.91% accordingly. This is consistent with the general fact that ice accretion will decrease the flight performance of the aircraft.

The estimated stability boundary provides limitation of the aircraft state parameters under icing conditions. When conducting envelope protection, the flight control should ensure the flight state within the area enclosed by the estimated stability boundary. If the aircraft is perturbed by the outside disturbance, as long as the flight state is within this boundary, the aircraft will recover to its trim point. Once the flight state exceeds the limitation due to some adverse conditions, for example, wind gust, the pilot or flight control should take action to prevent the aircraft entering into an irreversible dangerous situation.

##### 5.2. Validation of the ROA Results

To verify the correctness of the calculated ROA boundary, we use numerical simulation methods to judge if the method is acceptable or reliable. As stated above, if every flight state inside the calculated stability boundary can converge to the trim condition and area of the ROA is large enough to contain most of the stability area, the calculated ROA boundary can be considered effective.

In this paper, we use nonlinear conventional short period model in (9) to simulate the response of the aircraft starting from many initial state with inputs held fixed at their trim value . Projection of the simulation results onto* α*-

*q*plane along with the estimated stability boundary curve are shown in Figure 5(a). Red lines represent the unstable trajectories, blue lines represent the stable trajectories, and the bold green curve denotes the estimated stability boundary, that is, the ROA bound. It can be obviously seen that all trajectories inside the calculated ROA are stable and none of the unstable trajectory enter into the ROA.

**(a) clean**

**(b)**

**(c)**

**(d)**Each estimated stability boundary of these icing conditions is also presented on* α*-

*q*phase plane of the conventional short period model as shown in Figures 5(b), 5(c), and 5(d), respectively. As can be seen in Figure 5, whatever the icing severity is, the calculated ROA occupies most of the stability area of conventional short period model in the valid range of polynomial model. Therefore, it is feasible of using polynomial model to estimate the iced ROA of the aircraft, and the boundary of the ROA in different icing conditions can be approximately considered as the stability boundary of the aircraft.

According to the above approximation procedure, accuracy of the iced aerodynamic coefficients has great influence on the ultimate results of the stability boundary. The established nonlinear icing effect model in this paper, which is based on the existing wind tunnel data of a typical aircraft, can reflect the aerodynamic variation trend for a wide range of angle of attack under icing conditions. It is a proper model to verify the effectiveness of computing method of the iced stability boundary and can be used to research the iced nonlinear flight dynamic characteristics. However, precision iced aerodynamic coefficients for a specific aircraft must be available if a reliable stability boundary is needed for onboard ice protection system.

#### 6. Conclusions

A method has been developed for constructing the stability boundary of the aircraft under icing conditions through analyzing the region of attraction around the equilibrium point. Nonlinear icing effect model and polynomial model of the longitudinal dynamic system are constructed to estimate the ROA under icing conditions. The calculated results show that ice will deteriorate the stability characteristics of the aircraft, and the stability domain becomes less and less as icing severity increase.

Researches in this paper only involve two-dimensional ROA analysis and stability boundary estimation. Actually, as to three-dimensional or higher dimensional question, the method is also usable and the calculated stability boundary will be a multidimensional closed surface.

Ice accretion will not only cause degradation of the control effectiveness but also even disable the function of the control surface. This extreme condition, however, is not taken into consideration in this paper. To consider the control authority limit, further work should be done in the future.

#### Appendix

The polynomial period models and corresponding ROA in different icing conditions are presented as follows.

*(A1) *. The modified polynomial period model is

The calculated ROA is

*(A2) *. The modified polynomial period model is

The calculated ROA is

*(A3) *. The modified polynomial period model is

The calculated ROA is

#### Nomenclature

, , : | x and z body axes force and pitch moment coefficients |

L, D, M: | Lift force, drag force, and pitch moment, N |

U: | Equivalent air speed, m/s |

: | Angle of attack, rad |

q: | Pitch rate, body axis, rad/s |

: | Pitch angle, rad |

: | Normalized pitch rate, body axis, rad/s |

: | Mean aerodynamic chord, m |

: | Elevator deflection, rad |

: | Throttle position, 0~1. |

#### Competing Interests

The authors declare that they have no competing interests.

#### Acknowledgments

This study was cosupported by the State Key Development Program for Basic Research of China (no. 2015CB755800) and the National Natural Science Foundation of China (no. 61374145).