Aircraft nonlinear stability analysis and multidimensional stability region estimation under icing conditions

Liang QU,Yinghui LI,Haojun XU,Dengcheng ZHANG,Guoqiang YUAN

School of Aeronautics and Astronautics Engineering,Air Force Engineering University,Xi’an 710038,China

Aircraft nonlinear stability analysis and multidimensional stability region estimation under icing conditions

Liang QU,Yinghui LI*,Haojun XU,Dengcheng ZHANG,Guoqiang YUAN

School of Aeronautics and Astronautics Engineering,Air Force Engineering University,Xi’an 710038,China

Available online 14 February 2017

*Corresponding author.

E-mail address:liyinghui66@163.com(Y.LI).

Peer review under responsibility of Editorial Committee of CJA.

Production and hosting by Elsevier

http://dx.doi.org/10.1016/j.cja.2017.02.003

1000-9361©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.

This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Icing is one of the crucial factors that could pose great threat to flight safety,and thus research on stability and stability region of aircraft safety under icing conditions is significant for control and flight.Nonlinear dynamical equations and models of aerodynamic coefficients of an aircraft are set up in this paper to study the stability and stability region of the aircraft under an icing condition.Firstly,the equilibrium points of the iced aircraft system are calculated and analyzed based on the theory of differential equation stability.Secondly,according to the correlation theory about equilibrium points and the stability region,this paper estimates the multidimensional stability region of the aircraft,based on which the stability regions before and after icing are compared.Finally,the results are confirmed by the time history analysis.The results can give a reference for stability analysis and envelope protection of the nonlinear system of an iced aircraft.

©2017 Chinese Society of Aeronautics and Astronautics.Production and hosting by Elsevier Ltd.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

Dynamic system;

Equilibrium points;Icing aircraft;

Nonlinear stability;

Stability region

1.Introduction

Aircraft icing not only increases the weight of an aircraft,but also degrades both the performance and control of the aircraft because icing disrupts the flow of air over the aircraft1.Excessive icing accretion can lead to flow separation.It not only causes more resistance and less lift,but also leads to loss of control effectiveness2,or even aircraft crashes.Twelve percent of flight accidents were caused by icing according to the statistics about flight accidents caused by weather factors between 1999 and 2000,given by America Safety Advisor3.

On June 3rd 2006,a transport plane crashed in Anhui province.The direct cause is that the aircraft f l ew through ice clouds many times,and the freezing condition eventually drove the plane out of control.Incidents such as the American Eagle roll upset near Roselawn,Indiana in October 1994 and China Eastern Airlines CRJ-200 plane crash during takeoff at Baotou Inner Mongolia Airport in November 2004 are just two other accidents caused by aircraft icing4.The primary cause of these accidents is the effect of ice on the aircraft control effectiveness.

The effect of ice accretion on aircraft has been observed and investigated systematically since the 1940s5.In the early stage,the studies mainly focused on the changes of aerodynamic force and moment of an aircraft under icing conditions6–8.Afterwards,the performance,stability,and controllability of an aircraft after icing became the main interests9–11.Finally,the methodologies of envelope protection12,13,anti-icing14,and de-icing15of an aircraft under icing conditions were studied intensively.

The development of safer and more reliable aircraft must include better solutions for flights in icing and other bad weather conditions.Icing accidents can be prevented in two ways:(A)icing conditions can be avoided;(B)aircraft systems can be designed and operated in an ice-tolerant manner.For all kinds of aircraft,ice avoidance is the most desirable goal for enhancing safety and reliability.The main methods are:anti-icing coating,electro-thermal de-icing,electro-impulse de-icing,hot-air anti-icing,and so on.

However,for some aircraft,due to the funding or the weight constraint,ice tolerance will continue to be the preferred method.Sometimes,ice tolerance is necessary in some severe icing conditions.The research of this paper is on the aircraft stability and stability region under icing conditions,especially the multidimensional stability region.The stability region of nonlinear two-dimensional systems has been analyzed in some Refs.16,17.However,as we know,the twodimensional stability region solution is relatively simple and the stability of a two-parameter reaction aircraft system is not comprehensive.If a multidimensional stability boundary of the aircraft system under icing conditions can be obtained through analysis and calculation,we can know the aircraft parameters’tolerance in a certain degree of icing.Then,the aircraft’s safety can be judged clearly.

2.Model of aircraft dynamics

The purpose of this work is to analyze the equations of motions for an aircraft under an icing condition.The high angle of attack dynamics of the aircraft are inherently nonlinear.The nonlinearity will be stronger under the icing condition.Therefore,it is essential to study the stability performance of the iced aircraft.The stalling angle and available angle of attack which is the limit angle of attack for safety flight will decrease after icing as shown in Fig.1.Based on the above reasons,this work concentrates on the high angle of attack dynamics of an aircraft under an icing condition.

2.1.Equations of motion

The equations of motion used in this study assume a rigid aircraft.At the same time,the gyroscopic moment of the engine rotor is negligible,no applied thrust.The equations,consisting of translational equations and rotational equations18,19,are written in a principal axis system for analyzing the high angle of attack dynamics as follows:

wherev,g,andMare the aircraft velocity,gravity,and mass,respectively;θ is the pitch angle;α is the angle of attack;β is the angle of sideslip;φ is the roll angle;Sis the reference wing surface area;bis the reference wing span;¯cis the mean aerodynamic chord;p,q,andrare the body axis roll,pitch,and yaw rates,respectively;Ix,Iy,andIzare the moment inertias along aircraftxaxis,yaxis,andzaxis,respectively;Ixzis the product moment of inertia alongxandyaxes;Cxis the axial force coefficients;Cyis the lateral force coefficients;Czis the normal force coefficients;Clis the rolling moment coefficient;Cmis the pitching moment coefficients;Cnis the yawing moment coefficients;Q=0.5ρv2is the dynamic pressure,where ρ is the atmospheric density.

The change in aircraft velocity and gravity can be ignored since it has little impact on the equilibrium curve20when analyzing the high angle of attack dynamics.As a result,the tangential force equation can be discarded.Then the system can be reduced from a six-order system to a five-order system.

2.2.Aerodynamics force and moment

The aerodynamic force and moment models18are nonlinearities functions of multiple motion parameters,such as the angle of attack,the angle of sideslip,the elevator deflection,and so on.The main motion parameters are chosen to simplify the functions as follows:

where δeis the elevator deflection;δais the aileron deflection;δris the rudder deflection;Cyδa,Cyδr,Clδa,Clδr,Cnδa,andCnδrare control derivatives;Cxq,Cyp,Cyr,Czq,Clp,Clr,Cmq,Cnp,andCnrare stability derivatives.

One of the important characteristics of these aerodynamic parameters is that the data is discrete.To analyze aircraft motion characteristics with motion parameters’continuous transformation,the aerodynamic data is approximated with cubic interpolation in simulations.

2.3.Iced aircraft aerodynamics model

A simple,but effective representative icing model on aircraft flight mechanics was used in this paper,which was described in Ref.21.The icing severity factor can be represented by a linear variation parameter.The icing model is based on the following equation:

where η is the aircraft icing parameter,which represents the amount and severity of icing encountered in the given flying conditions of a particular aircraft for operation in the linear aerodynamic regime.KC（A）represents the change in an aircraft parameterC（A）,which is constant for a given aircraft.C（A）could represent any arbitrary performance,stability,or control derivative that is affected by ice accretion,andC（A）icedis its value under an icing condition.

It is assumed that the maximum possible value of the icing parameter η is denoted by ηmaxin a given condition.Hence, η can attain values between 0 and ηmaxin that condition.In other words,η is a varying parameter that increases from η =0 when the aircraft is clean,to a fully iced value of ηmax.A typical value of ηmaxis 1.The details on the variation of aircraft icing parameters can be referred to Ref.21.

3.Determination of equilibrium points and its analysis

The dynamics of the iced aircraft have been studied through determining the steady states of the equations of motion.The steady states are plotted as functions of the elevator deflection.

The main structural parameters of the example aircraft of this paper are as follows:M=9295.44 kg,S=27.87 m2,b=9.144 m,¯c=3.45 m.

The initial trim conditions are:altitudeH=4000 m,velocityv=70 m/s,aircraft icing parameter η=0.2.When fixing the aileron deflection δa=0°,the rudder deflection δr=-3°.Fig.2 shows the steady states curves of the aircraft with the elevator deflection under an icing condition,and the value ofKC（A）of the example aircraft can be referred to Ref.22.The stability of equilibrium points is analyzed by solving the Jacobi matrix of the system at the equilibrium points in the figure.

Fig.2 shows that when the elevator deflection is greater than-13°,the aircraft under an icing condition has stable equilibrium points(full line).It can be seen that the system has stable equilibrium points when the angle of attack is smaller than 25.8°in Fig.2(a),and at the same time,the modes of the iced aircraft are all stable.The system begins to show an unstable equilibrium solution since the angle of attack is beyond 25.8°.In this condition,it can be found that the Dutch roll model is unstable through calculating the characteristics roots of the system.

With the angle of attack continues to increase,the system enters another unstable state that the short period mode is unstable.We can see that the angle of attack maintains 38.7°in these unstable equilibrium points in Fig.2(a).The angle of attack is just the stall angle of attack of the aircraft in the icing condition.If the angle of attack continues to increase,the system will enter another complex unstable state.

The steady states of the angle of sideslip,the roll rate,the pitch rate,and the yaw rate shown in Fig.2(b)to(e)correspond to the steady state of the angle of attack shown in Fig.2(a).

Fig.3 shows the steady state of the angle of attack of the aircraft without icing in the same initial conditions.Comparing Fig.2(a)and Fig.3,we can find that the aircraft has a smaller stable angle of attack in the icing condition.The extent from a stable point to an unstable point in the same elevator deflection decreases.This illustrates that the stability region will decrease when the aircraft is under an icing condition.

4.Estimation of stability region

The stability analysis about the equilibrium point in Section 3 is based on the Lyapunov sense.Thus it’s only locally asymptotically stable.If the aircraft has many stable attractors,it is not enough to know the local asymptotic stability,because if the stability region of the expected stable equilibrium point is small,it is stable in theory,but its performance is unstable in reality.On the other hand,the stability analysis can’t indicate directly the aircraft’s stability margin or whether the aircraft will lose stability in a certain state.Consequently,it is very important to estimate the stability region of the aircraft’s stable equilibrium.The significance of the stability region is to characterize the anti-disturbance ability of the aircraft when at a stable state.The stability region decreases after icing,which means that the anti-disturbance ability is reduced,so it is necessary to give the multidimensional stability region.

Research about the nonlinear stability region of an aircraft is fairly limited.Refs.20,23used the bifurcation and catastrophe theory methodology to analyze nonlinear stability of an aircraft,and calculated the time history of the state parameters to know the aircraft’s motion.However,they didn’t study the stability region and stability boundary of the aircraft.

4.1.Method introduction

There are many methods available for estimating the stability region of a stable equilibrium point for nonlinear systems24,25.However,it is very difficult to apply these methods to estimate the stability region of a five-order system.Three parameters are picked out according to experience,which are most sensitive to the equilibrium points of a high angle of attack system.They are the angle of attack α,the angle of sideslip β,and the yaw rater.Hence,the stability region is defined in a threedimensional space of α,β,andr.

For a nonlinear autonomous system,

wherex（0）=0，x∈ Rn，f:Rn→ Rn.If~xmeetsf（~x）=0,~xis an equilibrium point of this system.LetJ（x）be the Jacobi matrix off（x）in~x,ifJ（x）has a positive real part characteristic root,~xis calledIunstable equilibrium point(UEP).IfJ（x）hasnpositive real part characteristic roots,~xis callednUEPs.IfJ（x）has no positive real part characteristic root,~xis called stable equilibrium point(SEP).IfJ（x）has no zero real part characteristic root,~xis called hyperbolic equilibrium point.

x（t）is the solution curve of Eq.(14),and for a hyperbolic stable equilibrium point~xof the system,its stable and unstable manifolds are:

Theorem 1.For system Eq.(14),~xis an UEP on the stability boundary ∂A（xs）of SEPxs,if and only ifwu（~x）∩A（xs）≠0,whereA（xs）is the stability region of the system24,26.

Theorem 2.Letxibe an UEP on the stability boundary of SEPxs,and then

whereEis the set of all SEPs24,26.

Theorem 1 indicates that if one of the unstable manifold of an UEP is beginning to become SEPxs,thenximust be on the stability boundary ofxs.Theorem 2 indicates that the stability boundary ofxsis the union of the stable manifold of UEPxion the stability boundary of.

The domain of attraction of the aircraft with a controller is obtained using the method of simulations in Refs.28,29.However,it is difficult to carry out simulations with every point in the state space because of the enormous computational database.According to Theorem 1,Theorem 2,and the waveform relaxation method of calculating the two-dimensional stability region of nonlinear dynamic systems30,the threedimensional stability region solution method can be summarized as follows:

(1)Find out all equilibrium points,analyze which point is an SEP,and determine which UEP is in the boundary of the asymptotic stability region of the stable equilibrium point.

(2)Calculate the two characteristics vectors u1and u2,which correspond to the two negative characteristic roots of the Jacobi matrix of the system in the UEP.

(3)In the plane consisted by u1and u2,a number of points are distributed uniformly on the circle with the UEP as its center and the smallerRas the radius of the initial points.Ris the radius of the initial circle for unstable manifold estimation.Use the negative system of nonlinear differential dynamic systems to f i x the stable manifolds of the UEP setting in the boundary.

The boundary of the asymptotic stability region of the SEP is composed of these stable manifolds.Fig.4 is the schematic diagram of the above method.

4.2.Multidimensional stability region of aircraft

The application of the method is described in the following iced aircraft system.

According to the results of the last section,the system of the iced aircraft has three equilibrium points when the elevator deflection is δe=-8°.Their states are shown in Table 1.

PointAis an SEP,while pointsBandCare UEPs.It can be ascertained that pointBis the point on the boundary of the asymptotic stability region of stable equilibrium pointA(it can be seen easily from Fig.2).Through the analysis of the characteristic roots of the Jacobi matrix of the nonlinear dynamic system in pointB,we can get that one of the three characteristic roots corresponding to α, β,andris positive,and the other two are negative.

Firstly,calculate the unit length stable eigenvectors corresponding to the two negative characteristic roots.Secondly,find out the intersection points of the spherical surface with the center at pointBand the radius equal to 0.1 and the surface composited with the two stable eigenvectors.Thirdly,calculate the integral trajectories of the negative system from these intersection points with the method in the last section.Finally,these trajectories are the stable manifolds of the stable equilibrium pointA.The curved surface consisting of these trajectories is the stability boundary.

The stability region obtained by the above method is described by the three-dimensional manifold around the stable equilibrium point,and in order to clearly distinguish the stability boundary,it is projected onto the twodimensional plane.Fig.5 shows the stability region of the stable equilibrium pointA.Each part of the figure shows the stability boundary in the α, β plane for different values ofr.It can be seen in the figure that these boundaries demonstrate clearly that a large region lies inside the stability region.Note that only part of the boundaries near the unstable equilibrium pointBis shown in the figure.The reason is that the manifold trajectories will be inaccurate with the increasing distance from pointBbecause of the sensitivity of the method to the accuracy of the initial point.However,the result can be a good judge of whether the aircraft will lose stability or not.

It can be seen that the stability region will become smaller and smaller with the yaw rate decreasing in Fig.5,which can be easily explained by the reason that the boundary is getting increasingly farther from the stable equilibrium points.

Fig.6 shows a comparison between the stability region boundaries in the α,β plane withr=-5.73（°）/s when the aircraft under an icing condition and a clean condition.It can be seen that the stability region will become smaller when the aircraft is under the icing condition,so the safety and stability of the iced aircraft need to be extra focused.

Table 1 States of equilibrium points(δe=-8°).

4.3.Time domain response verification

The time domain simulation from the initial point(α =35.5°,β =0°,r=-5.73（°）/s)in the stability region is shown in Fig.7.It can be seen that the state is stabilized at pointA.The time domain simulation from the initial point(α =35.5°,β =4°,r=-5.73（°）/s)out of the stability region is shown in Fig.8.It can be seen that the state diverges and is away from the stability point finally.These results verify the accuracy of the stability region given by Fig.5.

5.Conclusions

A nonlinear model of aircraft dynamics under icing conditions is established based on the equations of motion and aerodynamics force and moment affected by icing.The equilibrium surface of the system is obtained by solving the dynamic equations.Then the stability of equilibrium points is analyzed by solving the Jacobi matrix of the system at the equilibrium points.Finally,the three-dimensional stability region is obtained and verified by the time domain responses.

(1)The aircraft has a smaller stable angle of attack in an icing condition by solving the equilibrium points of the aircraft system before and after icing.Analyzing the state property of the aircraft yields the stability of the aircraft,but the anti-disturbance ability of the aircraft in a certain state cannot be obtained just by analyzing the stability.

(2)The multidimensional stability region of the nonlinear system of the aircraft is well estimated in the proposed method of this paper.Time domain simulation further verifies the accuracy of the presented method.The comparison results show that the stability region will become smaller when the aircraft is under an icing condition.

(3)The stability region clearly shows a visual parameter boundary and comprehensive information to pilots,so it is very useful to ensure aircraft safety and provide advice for operation when an aircraft flight encounters an icing condition.

Acknowledgements

This study was co-supported by the National Key Basic Research Program of China(No.2015CB755805)and the National Natural Science Foundation of China (No.61374145).

1.Chen Y,Luo YC,Liu GQ.Aircraft icing effect on flight safe.Equipment Manufact Technol2014;5:254–5[Chinese].

2.Reehorst A,Chung J,Potapczuk M,Choo Y,Wright W.An experimental and numerical study of icing effects on the performance and controllability of a twin engine aircraft.37th AIAA aerospace sciences meetingamp;exhibit.Reston:AIAA;1999.

3.Safety Advisor.Aircraft icing[Internet].2013 May[cited 2016 April 15].Available from:http://www.aopa.org/-/media/Files/AOPA/Home/Pilot%20Resources/ASI/Safety%20Advisors/sall.pdf.

4.Whalen E,Bragg MB.Aircraft characterization in icing using flight-test data.J Aircraft2005;42(3):792–4.

5.Kanter M.Flight performance on XB-25E aircraft in natural ice during February,March and April 1945.Washington,D.C.:Army Air Force;1945,Report No.:42-32281.

6.Gray VH,Von Glahn UH.Effect of ice and frost formations on drag of NACA 65(1)-212 airfoil for various modes of thermal ice protection.Washington,D.C.:NACA;1953,Report No.:NACATN-2962..

7.Gray VH,Von Glahn UH.Aerodynamic effects caused by icing of an unswept NACA 65A004 airfoil.Washington,D.C.:NACA;1958,Report No.:NACA-TN-4155.

8.Yi X,Gui YW,Zhu GL.Numerical method of a threedimensional ice accretion model of aircraft.Acta Aeronaut Astronaut Sinica2010;31(11):2152–8[Chinese].

9.Ranaudo RJ.Performance degradation of a typical twin engine commuter type aircraft in measured natural icing conditionsAIAA 22nd aerospace sciences meeting.Reston:AIAA;1984.

10.Wang MF,Wang LX,Huang CT.Computational effects of ice accretion on aircraft longitudinal stability and control.J Beijing Univ Aeronaut Astronaut2008;34(5):592–5[Chinese].

11.Zhang Q,Liu Y,Gao ZH.Simulation of aircraft flight dynamics affected by ice accretion.Flight Dynam2011;29(3):4–7[Chinese].

12.Bragg MB,Basar T,Perkins WR,Selig M,Voulgaris P.Smart icing systems for aircraft icing safety.Reston:AIAA;2002,Report No.:AIAA-2002–0813.

13.Sharma V,Voulgaris PG,Frazzoli E.Aircraft autopilot analysis and envelope protection for operation under icing conditions.J Guid Control Dynam2004;27(3):454–65.

14.Xiao CH,Zhang W,Gui YW.Test study on anti-icing effects of hydrophobic coating in icing wind tunnel.J Exp Fluid Mech2013;27(2):41–5[Chinese].

15.Li GC,He J,Lin GP.Electro-impulse de-icing(EIDI)technology study.J Aerospace Power2011;26(8):1728–35[Chinese].

16.Fujisaki Y,Sakuwa R.Estimation of asymptotic stability regions via homogeneous polynomial Lyapunov functions.SICE Annual Conference in Fukui.Piscataway(NJ):IEEE Press;2003.p.1411–4.

17.Zou Y,Yao H.Stable field analysis of aircraft flying at high angle of attack.J Air Force Eng Univ(Nat Sci Ed)2006;7(4):11–3[Chinese].

18.Jahnke CC,Culick FEC.Application of bifurcation theory to the high-angle-of-attack dynamics of the F-14.J Aircraft1994;31(1):26–34.

19.Gao H,Zhou ZQ.A study of the global stability of high performance aircrafts at high angle-of-attack.Acta Aeronaut Astronaut Sinica1987;8(11):562–71[Chinese].

20.He ZD,Guo W.A novel synthetical method for study nonlinear flight stability.Acta Aerodynam Sinica1990;8(2):143–51[Chinese].

21.Bragg MB,Hutchison T,Merret J.Effect of ice accretion on aircraft flight dynamics38th AIAA aerospace sciences meetingamp;exhibit.Reston:AIAA;2000.

22.Lampton A,Valasek J.Prediction of icing effects on the coupled dynamic response of light airplanes.J Guid Control Dynam2008;31(3):656–73.

23.Li K,Fang ZP.Application of bifurcation analysis to control law design at high angles of attack.Acta Aeronaut Astronaut Sinica2003;24(4):289–92[Chinese].

24.Vannelli A,Vidyasagar M.Maximal Lyapunov functions and domains of attraction for autonomous nonlinear systems.Automatica1985;21(1):69–80.

25.Chiang H,Hirsch MW,Wu F.Stability regions of nonlinear autonomous dynamical systems.IEEE Trans Autom Control1988;33(1):16–27.

26.Chiang HD,Wu FF,Varoiya P.Foundation of direct method for power system transient stability analysis.IEEE Trans Circ Syst1987;34(2):712–28.

27.Qu L,Li YH,Xu HJ.Analysis of longitudinal nonlinear stabilizing region for icing aircraft based on phase-plane method.Acta Aeronaut Astronaut Sinica2016;37(3):865–72[Chinese].

28.Li K,Fang ZP.High angle-of-attack control law design based on global stability analysis.J Beijing Univ Aeronaut Astronaut2004;30(6):516–9[Chinese].

29.Amitabh S,Girish D.Synthesis of nonlinear controller to recover an unstable aircraft from poststall regime.J Guid Control Dynam1999;22(5):710–7.

30.Lin XL.The waveform relaxation method for estimating the asymptotic stability region of nonlinear differential dynamic systems.Chinese J Eng Math2010;27(3):479–86[Chinese].

15 April 2016;revised 26 August 2016;accepted 30 November 2016