# Lyapunov Based Estimation of Flight Stability Boundary under Icing Conditions

Authors: Binbin Pei; Haojun Xu; Yuan Xue
Journal: Mathematical Problems in Engineering
 (2017)

---

## 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.


---

## Body

## 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 H∞ 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 H∞ 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. Section2 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,x,z, and the pitch moment coefficients are calculated by(1)CZ=CZ0+CZαα+CZδeδe+CZqqc-2V,CX=CX0+CXαα+CXα2α2,Cm=Cm0+Cmαα+Cmδeδe+Cmqqc-2V.Lift coefficient and drag coefficient can be available via a rotation(2)CDCL=cosαsinα-sinαcosαCXCZ.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:(3)CAiced=1+ηicekCA′CA,where CA is any arbitrary aerodynamic coefficient that is affected by ice accretion. ηice represents icing severity which is only the function of the atmospheric conditions. kCA′ 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](4)CX=CXα+ΔCXα,δe+ΔCXα,q,CZ=CZα+ΔCZα,δe+ΔCZα,q,Cm=Cmα+ΔCmα,δe+ΔCmα,q,where C∗α, ΔC∗α,δe, and ΔC∗α,q 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 forCZ, the CZ0+CZαα term in linear aerodynamic model (see (1)) is corresponding to the CZα term in nonlinear aerodynamic model (see (1)). And CZδe, CZq terms are corresponding to ΔCZα,δe,ΔCZα,q terms, respectively. According to linear icing effect model shown as (3), each term of CZ in the nonlinear aerodynamics model can be calculated by(5)CZαiced=CZαclean-CZ0clean1+kCzαη+CZ0clean1+kCz0η,ΔCZα,δeiced=CZα,δeclean1+kCzδeη,ΔCZα,qiced=CZα,qclean1+kCzqη,where 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 kCA. In fact, kCA 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,kCA is the same as linear icing effect model.(2)
In the stall region, value ofkCA 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 ofkCA 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, kCA turn out to be zero, which means the differences between clean and iced configuration disappear through point.Figure 1
Comparison between clean and iced aerodynamic coefficient with and without modification.
(a)
C L (b)
C D (c)
C mAs to drag coefficient termCDα, 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 CDα should be counted as separate problem. The computing method is also based on (3) but the value of kCD depends on the angle of attack.(6)CDiced=1+ηicekCDαCD.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α<αstall. However, when calculating the flight dynamic characteristics in stall and poststall area, the method is not proper. To solve this problem, the reduction factor εr is defined to measure the alteration of the stall angle. By the definition of εr, the original incidence angle interval can be mapped into a new interval as follows:(7)αCL′=εrαr,αr∈α0,αp,where εr∈0,1 and it has a negative correlation with icing severity. α0 is the initial point in the look-up table and αp 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.(8)αCm′=εeαe,αe∈α0,αp, where εe>1 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 Figure1. For brevity, only static parameters CLα,CDα,Cmα are shown.In Figure1, 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.


## 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,x,z, and the pitch moment coefficients are calculated by(1)CZ=CZ0+CZαα+CZδeδe+CZqqc-2V,CX=CX0+CXαα+CXα2α2,Cm=Cm0+Cmαα+Cmδeδe+Cmqqc-2V.Lift coefficient and drag coefficient can be available via a rotation(2)CDCL=cosαsinα-sinαcosαCXCZ.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:(3)CAiced=1+ηicekCA′CA,where CA is any arbitrary aerodynamic coefficient that is affected by ice accretion. ηice represents icing severity which is only the function of the atmospheric conditions. kCA′ 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](4)CX=CXα+ΔCXα,δe+ΔCXα,q,CZ=CZα+ΔCZα,δe+ΔCZα,q,Cm=Cmα+ΔCmα,δe+ΔCmα,q,where C∗α, ΔC∗α,δe, and ΔC∗α,q 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 forCZ, the CZ0+CZαα term in linear aerodynamic model (see (1)) is corresponding to the CZα term in nonlinear aerodynamic model (see (1)). And CZδe, CZq terms are corresponding to ΔCZα,δe,ΔCZα,q terms, respectively. According to linear icing effect model shown as (3), each term of CZ in the nonlinear aerodynamics model can be calculated by(5)CZαiced=CZαclean-CZ0clean1+kCzαη+CZ0clean1+kCz0η,ΔCZα,δeiced=CZα,δeclean1+kCzδeη,ΔCZα,qiced=CZα,qclean1+kCzqη,where 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 kCA. In fact, kCA 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,kCA is the same as linear icing effect model.(2)
In the stall region, value ofkCA 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 ofkCA 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, kCA turn out to be zero, which means the differences between clean and iced configuration disappear through point.Figure 1
Comparison between clean and iced aerodynamic coefficient with and without modification.
(a)
C L (b)
C D (c)
C mAs to drag coefficient termCDα, 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 CDα should be counted as separate problem. The computing method is also based on (3) but the value of kCD depends on the angle of attack.(6)CDiced=1+ηicekCDαCD.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α<αstall. However, when calculating the flight dynamic characteristics in stall and poststall area, the method is not proper. To solve this problem, the reduction factor εr is defined to measure the alteration of the stall angle. By the definition of εr, the original incidence angle interval can be mapped into a new interval as follows:(7)αCL′=εrαr,αr∈α0,αp,where εr∈0,1 and it has a negative correlation with icing severity. α0 is the initial point in the look-up table and αp 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.(8)αCm′=εeαe,αe∈α0,αp, where εe>1 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 Figure1. For brevity, only static parameters CLα,CDα,Cmα are shown.In Figure1, 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.


## 2.2.1. Establishment of the Model

The adopted equation form of the longitude nonlinear aerodynamics model is [2, 13](4)CX=CXα+ΔCXα,δe+ΔCXα,q,CZ=CZα+ΔCZα,δe+ΔCZα,q,Cm=Cmα+ΔCmα,δe+ΔCmα,q,where C∗α, ΔC∗α,δe, and ΔC∗α,q 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 forCZ, the CZ0+CZαα term in linear aerodynamic model (see (1)) is corresponding to the CZα term in nonlinear aerodynamic model (see (1)). And CZδe, CZq terms are corresponding to ΔCZα,δe,ΔCZα,q terms, respectively. According to linear icing effect model shown as (3), each term of CZ in the nonlinear aerodynamics model can be calculated by(5)CZαiced=CZαclean-CZ0clean1+kCzαη+CZ0clean1+kCz0η,ΔCZα,δeiced=CZα,δeclean1+kCzδeη,ΔCZα,qiced=CZα,qclean1+kCzqη,where 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 kCA. In fact, kCA 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,kCA is the same as linear icing effect model.(2)
In the stall region, value ofkCA 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 ofkCA 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, kCA turn out to be zero, which means the differences between clean and iced configuration disappear through point.Figure 1
Comparison between clean and iced aerodynamic coefficient with and without modification.
(a)
C L (b)
C D (c)
C mAs to drag coefficient termCDα, 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 CDα should be counted as separate problem. The computing method is also based on (3) but the value of kCD depends on the angle of attack.(6)CDiced=1+ηicekCDαCD.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α<αstall. However, when calculating the flight dynamic characteristics in stall and poststall area, the method is not proper. To solve this problem, the reduction factor εr is defined to measure the alteration of the stall angle. By the definition of εr, the original incidence angle interval can be mapped into a new interval as follows:(7)αCL′=εrαr,αr∈α0,αp,where εr∈0,1 and it has a negative correlation with icing severity. α0 is the initial point in the look-up table and αp 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.(8)αCm′=εeαe,αe∈α0,αp, where εe>1 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 Figure1. For brevity, only static parameters CLα,CDα,Cmα are shown.In Figure1, 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]:(9)U˙=1m-D-mgsin⁡θ-α+Txcos⁡α+Tzsin⁡α,α˙=1mU-L+mgcos⁡θ-α-Txsin⁡α+Tzcos⁡α+q,q˙=M+TmIyy,θ˙=q,wherem is mass,U is the airspeed, α is the angle of attack, θ is the pitch angle, and q is the pitch rate. Tx and Tz are projection of the thrust along the body axes x and z, respectively, and Tm denotes the pitch moment produced by thrust which is calculated by function of Tx and Tz. The lift forceL, drag forceD, and pitch momentM in the model are given by(10)L=q-SCLα,δe,q^,D=q-SCDα,δe,q^,M=q-Sc-Cmα,δe,q^,where q-≔1/2ρV2 is the dynamic pressure and q^=c-/2Vq is the normalized pitch rate (unitless). C∗ is unitless aerodynamic coefficient and c- 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:(11)sin⁡x≈x-16x3,cos⁡x≈1-12x2.Error is in a reasonable range forx≤50° 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.(12)U-1=aU2+bU+c,wherea,b, andc are constant coefficients to be determined.It is assume that the engine thrust is only relevant to the throttle positionδth; 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α andq) and control variables (like δe). 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 CLα, CDα, and CLα,q- can be given by(13)CLα=c0+c1α+c2α2+c3α3+c4α4,CDα=c0+c1α+c2α2+c3α3+c4α4+c5α5+c6α6,CLα,q-=k0+k1α+k2q-+k3q-2+k4q-α+k5α2,where ki,ci 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:(14)x˙=fx,u,where x≔V,α,q,θ and u∶=δe,δth.


#### 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 atH= 500 m andV= 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
Trim angle of attack and inputs versus icing severity.Figure2 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α,q are considered for simplicity.For polynomial model, the short period model can be extracted from the 4-state model by holdingV andθ and δth at their trim values Vt, θt, and δth,t:(15)α˙=f2Vt,α,q,θt,δth,t,δe,q˙=f3Vt,α,q,θt,δth,t,δe,where f2,f3 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(16)α˙=-2.1251α4+0.5278α3+4.1502α2+0.2427αδe+0.0027αq-0.0656δe2+0.0046q2-3.5905α-0.3254δe+0.9328q+0.1481,q˙=31.7299α3+0.8212q3-10.2418α2+48.8657αδe-59.9981α-73.0342δe-4.4190q+6.8423.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α0,q0 with the elevator inputs held fixed at its trim value. For comparison, trajectories of parametersα andq are shown on α-q phase plane. As the valid range of angle of attack is larger than −5°, only -5°≤α≤40° situations are demonstrated here as shown in Figure 3.Figure 3
Phase-plane comparison of polynomial short period and conventional short period models.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α andq 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.


## 3.1. Conventional Longitudinal Dynamic Model

The conventional longitudinal nonlinear dynamic model of an aircraft can be described as follows [18]:(9)U˙=1m-D-mgsin⁡θ-α+Txcos⁡α+Tzsin⁡α,α˙=1mU-L+mgcos⁡θ-α-Txsin⁡α+Tzcos⁡α+q,q˙=M+TmIyy,θ˙=q,wherem is mass,U is the airspeed, α is the angle of attack, θ is the pitch angle, and q is the pitch rate. Tx and Tz are projection of the thrust along the body axes x and z, respectively, and Tm denotes the pitch moment produced by thrust which is calculated by function of Tx and Tz. The lift forceL, drag forceD, and pitch momentM in the model are given by(10)L=q-SCLα,δe,q^,D=q-SCDα,δe,q^,M=q-Sc-Cmα,δe,q^,where q-≔1/2ρV2 is the dynamic pressure and q^=c-/2Vq is the normalized pitch rate (unitless). C∗ is unitless aerodynamic coefficient and c- 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:(11)sin⁡x≈x-16x3,cos⁡x≈1-12x2.Error is in a reasonable range forx≤50° 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.(12)U-1=aU2+bU+c,wherea,b, andc are constant coefficients to be determined.It is assume that the engine thrust is only relevant to the throttle positionδth; 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α andq) and control variables (like δe). 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 CLα, CDα, and CLα,q- can be given by(13)CLα=c0+c1α+c2α2+c3α3+c4α4,CDα=c0+c1α+c2α2+c3α3+c4α4+c5α5+c6α6,CLα,q-=k0+k1α+k2q-+k3q-2+k4q-α+k5α2,where ki,ci 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:(14)x˙=fx,u,where x≔V,α,q,θ and u∶=δe,δth.


### 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 atH= 500 m andV= 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
Trim angle of attack and inputs versus icing severity.Figure2 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α,q are considered for simplicity.For polynomial model, the short period model can be extracted from the 4-state model by holdingV andθ and δth at their trim values Vt, θt, and δth,t:(15)α˙=f2Vt,α,q,θt,δth,t,δe,q˙=f3Vt,α,q,θt,δth,t,δe,where f2,f3 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(16)α˙=-2.1251α4+0.5278α3+4.1502α2+0.2427αδe+0.0027αq-0.0656δe2+0.0046q2-3.5905α-0.3254δe+0.9328q+0.1481,q˙=31.7299α3+0.8212q3-10.2418α2+48.8657αδe-59.9981α-73.0342δe-4.4190q+6.8423.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α0,q0 with the elevator inputs held fixed at its trim value. For comparison, trajectories of parametersα andq are shown on α-q phase plane. As the valid range of angle of attack is larger than −5°, only -5°≤α≤40° situations are demonstrated here as shown in Figure 3.Figure 3
Phase-plane comparison of polynomial short period and conventional short period models.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α andq 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.


## 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:(11)sin⁡x≈x-16x3,cos⁡x≈1-12x2.Error is in a reasonable range forx≤50° 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.(12)U-1=aU2+bU+c,wherea,b, andc are constant coefficients to be determined.It is assume that the engine thrust is only relevant to the throttle positionδth; 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α andq) and control variables (like δe). 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 CLα, CDα, and CLα,q- can be given by(13)CLα=c0+c1α+c2α2+c3α3+c4α4,CDα=c0+c1α+c2α2+c3α3+c4α4+c5α5+c6α6,CLα,q-=k0+k1α+k2q-+k3q-2+k4q-α+k5α2,where ki,ci 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:(14)x˙=fx,u,where x≔V,α,q,θ and u∶=δe,δth.


## 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 atH= 500 m andV= 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
Trim angle of attack and inputs versus icing severity.Figure2 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α,q are considered for simplicity.For polynomial model, the short period model can be extracted from the 4-state model by holdingV andθ and δth at their trim values Vt, θt, and δth,t:(15)α˙=f2Vt,α,q,θt,δth,t,δe,q˙=f3Vt,α,q,θt,δth,t,δe,where f2,f3 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(16)α˙=-2.1251α4+0.5278α3+4.1502α2+0.2427αδe+0.0027αq-0.0656δe2+0.0046q2-3.5905α-0.3254δe+0.9328q+0.1481,q˙=31.7299α3+0.8212q3-10.2418α2+48.8657αδe-59.9981α-73.0342δe-4.4190q+6.8423.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α0,q0 with the elevator inputs held fixed at its trim value. For comparison, trajectories of parametersα andq are shown on α-q phase plane. As the valid range of angle of attack is larger than −5°, only -5°≤α≤40° situations are demonstrated here as shown in Figure 3.Figure 3
Phase-plane comparison of polynomial short period and conventional short period models.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α andq 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](17)x˙=fx,x0=x0,where xt∈Rn is the state vector and f:Rn→Rn is a vector polynomial function of x with f0=0. Assume that the origin is a locally asymptotically stable equilibrium point; the ROA of the origin is defined as(18)Ω≔x0∈Rn∣limt→∞⁡φt,x0=0,where φt,x0 is the solution of (1) that starts at initial state x0.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γ>0 and a polynomial V∶Rn→R such that [18](1)
V0=0 and ∀x≠0,Vx>0;(2)
Ωγ≔x∈Rn∣Vx≤γ is bounded;(3)
Ωγ⊂x∈Rn∣∇Vxfx<0∪0,then for all x0∈Ωγ, the solution of (17) exists and limt→∞xt∈0. As such, Ωγ is a subset of the ROA for (17). A functionV 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:(19)εβ≔x∈Rn∣px≤βand it maximizes εβ under the constraint εβ⊆Ωγ. β is a positive value and px is shape function which is normally defined as px∶=xTNx. The choice of px 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:(20)maxβoverV∈Pn,V0=0,l1,l2,s1,s2∈Σns.t.V-l1∈Σn-β-ps1+γ-V∈Σn-γ-Vs2-∂V∂xf+l2∈Σn,where Pn is the set of all polynomials inn variables and Σn is the set of SOS polynomials in n variables. l1,l2 are positive definite polynomial of the form(21)li=∑j=1nεijxj2for i=1,2 and εij are small positive constants on the order of 10-6. In this paper, V-s iteration algorithm [26] is adopted to solve the SOSP. Finally, the stability boundary ES of the polynomial nonlinear system can be depicted as (22)ES≔x∈Rn∣Vx=γ.
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 Section2. 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:(23)α˙=-2.1251α4+0.2127α3+4.1913α2+0.0027αq+0.0046q2-3.2655α+0.9329q,q˙=31.7299α3+0.8212q3-6.7138α2-57.4663α-4.4190q.Based on SOS theory, ROA around the origin point is(24)Ωγ≔α,q∈R∣0.8705α2-0.04338αq+0.0239q2≤0.1472.Thus the estimated stability boundary is(25)ES≔α,q∈R∣0.8705α2-0.04338αq+0.0239q2=0.1472.Similarly, the modified polynomial short period model and the corresponding ROA in different icing conditionsη=0.1, 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 -5°≤α≤50° due to the limitation of the data available. Thus, the effective region of the ROA is enclosed by vertical line α=-5° (shown in green dashed line) and right part of ES 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.Figure 4
Comparison of the stability boundary for different icing conditions.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 × 103, for iced configuration of η=0.1, 0.2, and 0.3, and the effective area of the ROA came out to be 5.5576 × 103, 3.6134 × 103, and 2.1872 × 103, 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 Vt,α0,q0,θt with inputs held fixed at their trim value δth,t,δe,t. 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.Figure 5
Phase plane ofα-q for clean and iced configuration of conventional short period model.
(a)
clean (b)
η = 0.1 (c)
η = 0.2 (d)
η = 0.3Each 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.


## 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 Section2. 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:(23)α˙=-2.1251α4+0.2127α3+4.1913α2+0.0027αq+0.0046q2-3.2655α+0.9329q,q˙=31.7299α3+0.8212q3-6.7138α2-57.4663α-4.4190q.Based on SOS theory, ROA around the origin point is(24)Ωγ≔α,q∈R∣0.8705α2-0.04338αq+0.0239q2≤0.1472.Thus the estimated stability boundary is(25)ES≔α,q∈R∣0.8705α2-0.04338αq+0.0239q2=0.1472.Similarly, the modified polynomial short period model and the corresponding ROA in different icing conditionsη=0.1, 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 -5°≤α≤50° due to the limitation of the data available. Thus, the effective region of the ROA is enclosed by vertical line α=-5° (shown in green dashed line) and right part of ES 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.Figure 4
Comparison of the stability boundary for different icing conditions.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 × 103, for iced configuration of η=0.1, 0.2, and 0.3, and the effective area of the ROA came out to be 5.5576 × 103, 3.6134 × 103, and 2.1872 × 103, 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 Vt,α0,q0,θt with inputs held fixed at their trim value δth,t,δe,t. 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.Figure 5
Phase plane ofα-q for clean and iced configuration of conventional short period model.
(a)
clean (b)
η = 0.1 (c)
η = 0.2 (d)
η = 0.3Each 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.



---

*Source: 6901894-2017-01-30.xml*