Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 1, pp. 269-284, Warsaw 2012 50th anniversary of JTAM AIRCRAFT DYNAMICS DURING FLIGHT IN ICING CONDITIONS Grzegorz Kowaleczko Michał Wachłaczenko Air Force Institute of Technology, Warsaw, Poland e-mail: g.kowaleczko@chello.pl; michal.wachlaczenko@itwl.pl This paper presents numerical analysis of airplane motion parameters during wing in-flight icing. The wing external flaps deflection event was taken under consideration.Themathematicalmodel of the airplanemo- tion including twelveordinarydifferential equationswas employed.Using available publications as well as TS-11 “Iskra” jet trainer geometrical and mass relationships, a computer code evaluating motion parameters of an iced airplanewas developed. On the basis of analysis results, conc- lusions referring airplane dynamic properties during wing ice accretion was formulated. Key words: ice-accretion, wing icing, simulation Notations α,β – angle of attack and slideslip Φ,Θ,Ψ – aircraft roll, pitch and yaw angle m – airplane total mass V = [u,v,w] ⊤ – centre of mass of airplane velocity vector Ω= [p,q,r] ⊤ – rotational velocity vector in airplane-fixed system of coordi- nates F = [X,Y,Z] ⊤ – external forces vector acting on the airplane M = [L,M,N] ⊤ – moments of forces vector M acting on the airplane K = ∑ i (ri×miV i) – total angular momentum around taken reference point Mgir = Jω×Ω – gyrostatic moment 270 G. Kowaleczko, M. Wachłaczenko ω= [ω,0,0] ⊤ – engine rotor angular velocity vector J – engine rotor moment of inertia CDa,Cya,CLa – aircraft drag, side and lift coefficients Cl,Cm,Cn – aircraft roll, pitch, yaw moment coefficients bA – mean wing chord ρ – air density l,S – wing span and wing area t0, t – beginning of icing phenomenon and flight current time Cobla (α,t) – aerodynamic coefficient current value of iced airplane Ca(α,t0) – aerodynamic coefficient value of non-iced airplane ∆Ca(α,t) – change of aerodynamic coefficient value caused by icing 1. Introduction The airplane in-flight icing is a serious problem, still causing many accidents. It can proceed in every climatic zone in foster outer conditions. Usually, the changes of aerodynamic characteristics are so big that the aircrews cannot counteract the icing effects even while using highly efficient onboard anti- and de-icing equipment. Ice accretion on external airplane components is one of themost dangerous phenomenon which can take place during flight in variable earth atmospheric conditions. It is known since the very beginning of world aviation, but still not fully investigated. Formation of a solid ice cover with different structures and shapes of different speed rates and intensities on airplane external surfa- ces causes an increase of roughness and overall airplane mass. The process of deformation of an airfoil section of the wing and control surfaces by ice accre- tion has direct influence on dynamic properties, which cause flowdisturbances around the airplane, leading to variation of its dynamic characteristics. The risk, that is carried during flight in icing conditions, guided many re- search centers to gain theirs interest in this phenomenon. Conducted studies, theoretical and experimental, are aiming at as close as possible determination of ice accretion influence on aerodynamics and airplane dynamic properties. Investigationsmade in especially preparedwind tunnels andduring test-flights in icing conditions require major staff engagement and big expenses, so that onlywealthy research centers can afford it. Because of this, analytical conside- rations are made using specialized computer codes allowing one to determine Aircraft dynamics during flight in icing conditions 271 iced airfoils and airplane models aerodynamic characteristics. Nevertheless, thismethod of scientific research not always shows true airplane dynamic pro- perties during flight in icing conditions, that is why final conclusions must be based on an experimental research conducted on real object flights in weather conditions bringing ice accretion on airplane outer parts. Ice accretion can formonvarious airplane external elements.Depending on where firm ice cover appears (f.e. wing, tailplane, fuselage, etc.), the airplane may respond tomotion parameters changes and control displacements in diffe- rent ways as its aerodynamic balance is changing due to in-flight icing. In this research, it is assumed that icing leads to a decrease of aerodynamic lift, incre- ase of aerodynamic drag, pitching moment change, reduction of wing critical angle of attack value. The flap displacement produces a change in flow around horizontal tailplane. Asymmetrical wing icing cause airplane uncontrolled roll (asymmetric of aerodynamic lift) and yaw (asymmetric of aerodynamic drag). These assumptions were specified on the basis of own studies and literature data. Technical data of TS-11 “Iskra” jet trainer were obtained fromRewucki andSierputowski (2002) and aerodynamics characteristics were derived on the basis of Rley (1998). 2. Airplane equations of motion Spatialmotion of the airplane treated as a rigid body is described by a system of twelve nonlinear ordinary differential equations. In practice, engineering calculations employapproximatemethodsonthebasis of computer technology. Therearemanyexternal forces andmomentsof forces actingon theaircraft during flight, which are complex functions of motion parameters as well as geometrical and mass quantities of the investigated object. Additional forces andmoments of forces are generated by control surfaces displacements. The accepted simplifing assumptions: • airplane is a rigid bodywith constant mass, constant moment of inertia and invariable position of the centre of mass • control surfaces are stiff and their axes of rotation have an invariable position against the airplane • plane Oxz is a geometrical, mass and aerodynamic symmetry plane. The gyrostatic moment from revolving jet engine elements is considered in equations of motion. 272 G. Kowaleczko, M. Wachłaczenko Fig. 1. Systems of coordinates and angles defining its mutual position Themutual orientation of coordinate systems are defined by following an- gles: angle of attack α, angle of slideslip β, aircraft roll angle Φ, aircraft pitch angle Θ, aircraft yaw angle Ψ. Performing rotations successively by Euler angles Ψ,Θ and Φ, a Oxgygzg to Oxyz transition matrix is determined    x y z    =A    xg yg zg    where A=    cosΨ cosΘ sinΨ cosΘ −sinΘ cosΨ sinΘsinΦ− sinΨ cosΦ sinΨ sinΘsinΦ+cosΨ cosΦ cosΘsinΦ cosΨ sinΘcosΦ+sinΨ sinΦ sinΨ sinΘcosΦ− cosΨ sinΦ cosΘcosΦ    Performing rotations by −β and α angles, a Oxayaza to Oxyz transition matrix is determined    x y z    =    cosαcosβ −cosαsinβ −sinα sinβ cosβ 0 sinαcosβ −sinαsinβ cosα       xa ya za    2.1. Airplane translational equations of motion Vector equation of motion of the centre of mass m dV dt =F (2.1) Aircraft dynamics during flight in icing conditions 273 and it canbewritten inascalar formin theairplanefixedsystemof coordinates Oxyz m(u̇+qw−rv)=X m(v̇+ru−pw)=Y m(ẇ+pv−qu)=Z (2.2) Using transformations inKowaleczko (2003) andWachłaczenko (2010), the airplane translational equations of motion are obtained V̇ = 1 m (Xs cosβ+Ys sinβ) α̇= 1 cosβ− Zα̇ S mV [ZSst mV +qcosβ− (pcosα+r sinα)sinβ ] β̇= 1 mV (Yscosβ−Xs sinβ)+psinα−rcosα (2.3) Equation (2.2) considers that the aerodynamic force Pza as a part of for- ce Za depends, among other things, on the rate of change of the airplane angle of attack. So, it is calculated as follows Za =Zast+Z α̇ a α̇ (2.4) 2.2. Airplane rotary equations of motion The vector equation of the total angular momentum change dK dt =M +Mgir (2.5) Because the Oxz surface is the airplane geometrical, mass and aerodyna- mic symmetry plane and the engine rotor rotates with an angular velocity ω along axis parallel to airplane axis Ox, the scalar form of equation (2.5) in the airplane fixed system of coordinates Oxyz is Ixṗ− (Iy − Iz)qr− Ixz(ṙ+pq)=L Iyq̇− (Iz − Ix)pr− Ixz(r 2 −p2)=M−Jωr Izṙ− (Ix− Iy)pq− Ixz(ṗ−qr)=N+Jωq (2.6) 274 G. Kowaleczko, M. Wachłaczenko Converting system (2.6) (Kowaleczko, 2003) ṗ= 1 IXIZ − I 2 XZ { [L+(IY − IZ)qr+ IXZpq]IZ +[N +Jωq+(IX − IY )pq− IXZqr]IXZ } q̇= 1 IY [M−Jωr+(IZ − IX)rp+ IXZ(r 2 −p2)] ṙ= 1 IXIZ −I 2 XZ { [L+(IY − IZ)qr+IXZpq]IXZ +[N +Jωq+(IX − IY )pq− IXZqr]IX } (2.7) Moreover, to equations (2.3) and (2.7) we still need to add equations de- fining the current mass centre angular position (Euler angles, the kinematic relations) (2.8) and equations determining airplane centre of mass position components in gravitational system of coordinates Oxgygzg (2.9) Φ̇= p+(rcosΦ+q sinΦ)tanΘ Θ̇= qcosΦ−r sinΦ Ψ̇ =(rcosΦ+q sinΦ) 1 cosΘ (2.8) and ẋg =V [cosαcosβcosΘcosΨ+sinβ(sinΦsinΘcosΨ− cosΦsinΨ) +sinαcosβ(cosΦsinΘcosΨ+sinΦsinΨ)] ẏg =V [cosαcosβ cosΘsinΨ+sinβ(sinΦsinΘsinΨ+cosΦcosΨ) +sinαcosβ(cosΦsinΘsinΨ− sinΦcosΨ)] żg =V (−cosαcosβ sinΘ+sinβ sinΦcosΘ+sinαcosβcosΦcosΘ) (2.9) Equations (2.3), (2.7), (2.8) and (2.9) form a system of twelve nonlinear ordinary differential equations describing the airplane spatial motion conside- red as a rigid body with a constant mass. The motion parameters vector has the following components: V , α, β, p, q, r,Φ,Θ, Ψ, xg, yg, zg. 3. Forces and moments acting on the airplane The vector of external forces which are affecting the airplane (right-hand side of equation (2.1)) may be written as F =Q+T +R (3.1) Aircraft dynamics during flight in icing conditions 275 The airplane gravitational force Q in the Oxgygzg system has only one component Q= [0,0,mg]⊤. In Oxayaza system Qxa =mg(−cosαcosβ sinΘ+sinβ cosΘsinΦ+sinαcosβ cosΘcosΦ) Qya =mg(cosαsinβ sinΘ+cosβcosΘsinΦ− sinαsinβcosΘcosΦ) (3.2) Qza =mg(sinαsinΘ+cosαcosΘcosΦ) The engine thrust vector T lies in the airplane symmetry plane Oxz and it is parallel to Ox axis and placed in the airplane mass centre, so that in Oxyz system T = [T,0,0]⊤ (3.3) In the flow system of coordinates Txa =T cosαcosβ Tya =−T cosαsinβ Tza =−T sinα (3.4) The resultant aerodynamic R force projections on the Oxayaza system axes and the moment vector M = [L,M,N]⊤ on the right-hand side of equ- ation (2.5) is the resultant moment acting on the airplane Rxa =−Pxa =−CDa ρV 2 2 S L=Cl ρV 2 2 SbA Rya =−Pya =−Cya ρV 2 2 S M =Cm ρV 2 2 SbA Rza =−Pza =−CLa ρV 2 2 S N =Cn ρV 2 2 Sl (3.5) 4. Initial flight conditions At thebeginningof simulation, theairplaneperformsa rectilinear steadyflight with a given constant speed. For this condition, we can take p= q= r=0 β=0 Φ=0 γ =Θ−α=0 and V̇ = α̇= β̇= ṗ= q̇= ṙ= Θ̇= Φ̇=0 To perform a flight in these conditions, the resultant force and moment acting on the airplane must be equal to zero PDa sinα+(PLa−mg)cosα=0 T−PDacosα+(PLa−mg)sinα=0 (4.1) 276 G. Kowaleczko, M. Wachłaczenko Transformations of the above equations allow one to obtain a steady flight angle of attack α (4.1)1 and indispensable thrust T (4.1)2.Theknownairplane angle of attack αpermits one to set the elevator displacement angle δH0 (4.2), so as the pitching moment coefficient has to be equal to zero Cm =Cmst(V,α)+ ∂Cm ∂δH δH0 =0 δH0 =− Cmst ∂Cm ∂δH (4.2) 5. Aerodynamic force and moment coefficients The experimentally obtained airplane static aerodynamic characteristics of drag CDa(α,Ma), lift CLa(α,Ma)andpitch Cm(α,Ma)areused in this paper. Along with the determined dynamic derivatives, this allows one to set final expressions defining the force andmoment coefficients: – drag force: CDa =CDa(α,Ma) – side force: Cya =C β yaβ+C p yap+C r yar+C δV ya δV +C δl yaδl – lift force: CLa =CLa(α,Ma)+C α̇ Laα̇+C q La q+C δH La δH – rolling moment: Cl =C β l β+C p l p+Crl r+C δV l δV +C δl l δl – pitching moment: Cm =Cm(α,Ma)+C α̇ mα̇+C q mq+C δH m δH – yawing moment: Cn =C β nβ+C p np+C r nr+C δV n δV +C δl n δl 6. Airplane ice accretion effects modeling Nowadays, there aremanynumericalmodels employed to conduct studies abo- ut the aircraft surface ice accretion influence on aerodynamic characteristics and dynamic parameters. The computer code using equations of motion and their numerical integration allows one to study the object movement under the icing conditions in any given ice accretion models. The in-flight icing causes flow distortion around the lift and control surfa- ces, and produces considerable forces and moments of forces changes during flight. The change of airfoil geometry results in: • increase of aerodynamic drag (up to hundreds of percent) due to greater airplane surfaces roughness • decrease of the aerodynamic lift maximum value CLmax • decrease of the critical angle of attack value αcr Aircraft dynamics during flight in icing conditions 277 • decrease of the dCLa/dα derivative value • increase of the stall speed • additional pitching moment arising • asymmetric wing icing leads to uncontrolled airplane rolling and yawing. The assumption of this paper is that the relative value of the aerodynamic coefficient of forces andmoments of forces is changing as in Fig.2. The studies on the aircraft icing caused that in the numerical modeling, only the wing lift, drag and pitch coefficients changes are taken under consideration in the computer codes. Fig. 2. Aircraft maximum lift coefficient relative change as a function of time The aerodynamic coefficients change as follows Cobla (α,t) =Ca(α,t0)−∆Ca(α,t) (6.1) where ∆Ca denote the lift, drag and pitch coefficient, and they are functions of the airplane angle of attack and time (mathematical function elaborated by Wachłaczenko, beyond the scope of this paper). This model of airplane icing phenomenon does not include the total we- ight change (negligible increase, about 1.1% of clean airplane weight) and the displacement of mass centre (shift less then 0.01m) during flight in icing con- ditions. The moments of inertia and moments of deviation are constant as well. To study the influence of aircraft icing on the dynamic properties, using a computer code to determine flight parameters, different cases of aerodynamic coefficients change can be implemented.Themajority of scientific publications about airplane icing and its effect refer to decrease of critical angle of attack as well as aerodynamic lift, increase of total mass and aerodynamic drag. The 278 G. Kowaleczko, M. Wachłaczenko change of mass can be taken under account depending on the complexity of the used computer code determining iced airplane flight parameters, but the increase of total weight caused by ice accretion is only about few percent. The in-flight icing causes an aircraft balance change, hence the pitching moment coefficient alteration is taken under consideration. Iced airplane aerodynamic force andmoment of force coefficient measure- ments in especially prepared wind tunnels prove that in all cases the aerody- namic drag increases. It is caused by ice accretion on lifting surfaces, tailplane, fuselage and other airplane external equipment. The lift coefficient decreases depending on the type and shape of the ice cover on the wing, but it can also increase slightly. In the computer code which determines aircraft flight parameters, the ae- rodynamic coefficient change model is used (on the basis of equation (6.1)): – wing drag coefficient (Fig.3): CoblDa(α,t) =CDa(α,t0)+∆CDa(α,t) – wing lift coefficient (Fig.4): CoblLa(α,t) =CLa(α,t0)−∆CLa(α,t) – wing pitch coefficient (Fig.5): Coblm (α,t) =Cm(α,t0)−∆Cm(α,t) Themanner of pitching moment coefficient change is based on a very few experimental data, which the authors could find in the available literature. Figures 3-5 show the changes of airplane lift, drag, pitch coefficient altera- tions during ice-accretion simulations: ——————– no ice t< 10s —-—-—- – iced wing t=(10+30)s – – – – – – iced wing t=(10+60)s Fig. 3. Airplane drag coefficient Aircraft dynamics during flight in icing conditions 279 Fig. 4. Airplane lift coefficient Fig. 5. Airplane pitch coefficient 7. Aircraft in-flight icing analysis The initial condition is a steady flight. Depending on the calculation variant, it was a rectilinear flight (case a – solid line) and steady descent (case b – dot line) with velocity of 90m/s (324km/h, Ma ≈ 0.265) on altitude equal to 1000m. The airplane is performing flight in a smooth configuration (flaps on δKL = 0 ◦). The jet engine thrust was determined on the base of steady flight condition (equation (4.1)2) and remains constant in the analysed time interval. The static pitchingmoment is compensated by the horizontal control surface deflection (equation (4.2)). The system of twelve nonlinear ordinary differential equations (equations (2.3), (2.7)-(2.9)) and the wing icing model was used to perform the analysis. The calculations were carried for 200 seconds of flight. The in-flight icing phenomenon occurs in the 10th second of flight. There is no pilot’s response, so the airplane behaviour is investigated. 280 G. Kowaleczko, M. Wachłaczenko The symmetric wing icing starts in the 10th second of linear steady flight. The extra positive pitching moment (Fig.5) caused by ice accretion creates a slow increase in the airplane altitude (Fig.10). The airplane angle of attack (Fig.7) increases as well, and it approaches to its critical value which is lower than for an non-iced airplane. The temporary climb produces the loss of speed (Fig.6). After 50 seconds of flight, airplane has not got enough lift to continue the climb and begins to descent with a velocity increase. The angle of attack is maintained in the range of 8.0-8.5 degrees. Similar calculations were conducted for the initial descent case. In all the abovementioned events, the changes were alike. They only differ in themaxi- mum values of flight parameters. Figure 10 shows that in about two minutes from the start of the icing phenomenon, the airplane reaches zero altitude performing the right turn and deviating about 250 meters from the initial direction. Fig. 6. Airplane total velocity Fig. 7. Airplane angle of attack Aircraft dynamics during flight in icing conditions 281 Fig. 8. Airplane pitching angular velocity Fig. 9. Airplane angle of pitch Fig. 10. Airplane flight altitude 282 G. Kowaleczko, M. Wachłaczenko Fig. 11. Airplane trajectory projection on the ground plane 8. Conclusions The icing phenomenon is especially dangerous during takeoff and landing.Low velocity, high angles of attack close to its critical value, in such flight condi- tions the changes of aerodynamic lift generated by icing are the most consi- derable. Theconductedanalysis provide following conclusionsof aircrewperforming flights in icing conditions: • on the basis of meteorological information try to pass round the icing zone, in the case of entering – immediately leave icing zone • if there is a huge pitchingmoment while displacing the flaps, theymust be put in the previous position or fully hidden. If the in-flight icing on lifting or control surfaces occurs, the autopilot must be disabled (negligible changes of airplane dynamics would be leveled by autopilot and the aircrew would be unable to notice these changes as a result of possible icing). During the process of removing the ice cover from lifting surfaces by onboard deicing equipment, the process must be observed constantly, non-uniform tearing off of the ice cap can cause asymmetry of forces on the wing leading to uncontrolled airplane rolling and yawing. References 1. Al-Sharabi M., Maryniak J., 2002, Wyznaczanie charakterystyk aerody- namicznych oblodzonego samolotu, [In:] Mechanika w Lotnictwie, ML-X, J. Maryniak (Edit.), 121-134 Aircraft dynamics during flight in icing conditions 283 2. Björk Å., Dahlquist G., 1983,Metody numeryczne, PWN,Warszawa 3. Cebeci T., 1995, Effect of ice on airfoil stall at highReynolds numbers,AIAA Journal, 33, 7, 1351-1352 4. Dżygadło Z., Kowaleczko G., Sibilski K., 2001,Charakterystyki aerody- namiczne samolotu ze skrzydłem pasmowym i jego ruchy własne, Prace Insty- tutu Lotnictwa,Warszawa 5. Fiszdon W., 1962,Mechanika lotu Cz. I, PWN,Warszawa 6. Fiszdon W., 1962,Mechanika lotu Cz. II, PWN,Warszawa 7. FrantM.,KowaleczkoG., SobierajW., 2002,Doświadczalne charaktery- styki profilu lotniczego zmodelowanym oblodzeniem, Prace InstytutuLotnictwa, Warszawa 8. Gutowski R., 1971,Równania różniczkowe zwyczajne, WNT, Warszawa 9. Hersi I., Maryniak J., 1993, Symulacja numeryczna wpływu oblodzenia po- wierzchni nośnych samolotu na charakterystyki aerodynamiczne i parametry dy- namiczne, Prace Instytutu Lotnictwa,Warszawa 10. Kochański A., 1938,Obladzanie samolotów, Lwów 11. Kowaleczko G., 2003, Zagadnienie odwrotne w dynamice lotu statków po- wietrznych, WydawnictwoWAT,Warszawa 12. Kowaleczko G., Sobieraj W., 2002,Wpływ oblodzenia na charakterystyki lotne statków powietrznych, [In:]Mechanika w Lotnictwie, ML-X, J.Maryniak (Edit.), 87-112 13. Krzyżanowski A., 1984,Mechanika lotu, SkryptWAT,Warszawa 14. Maryniak J., Wojciechowski Z., 2002, Wpływ oblodzenia na charaktery- styki aerodynamicznewybranychmodeli samolotów, [In:]Mechanika w Lotnic- twie, ML-X, J. Maryniak (Edit.), 113-120 15. Ostosławskij N.W., 1957,Aerodinamika samoleta, Oborongizdat,Moskwa 16. PracaZbiorowa,1981,Wyniki badań tunelowych samolotuTS-11 „Iskra”,Opra- cowanie wewnętrzne ITLWAT,WAT,Warszawa 17. Rewucki P., Sierputowski P., 2002, Struktury oblodzenia skrzydła samo- lotu oraz ich wpływ na charakterystyki aerodynamiczne, [In:] Mechanika w Lotnictwie, ML-X, J. Maryniak (Edit.), 189-202 18. Riley J.T., 1998,Mixed-Phase Icing Conditions: AReview, Office ofAviation Research,Washington DC 19. Wachłaczenko M., 2010, The influence of wing transverse vibration on dy- namic parameters of an aircraft,Journal of Theoretical andAppliedMechanics, 48, 2, 397-413 284 G. Kowaleczko, M. Wachłaczenko Parametry lotu samolotu w warunkach oblodzenia Streszczenie Praca przedstawia wyniki numerycznej analizy parametrów ruchu samolotu pod- czas symetrycznego oblodzenia skrzydła oraz usterzenia wysokości w czasie lotu. Pod uwagę wzięto również przypadek wypuszczania klap zaskrzydłowychw czasie wystę- powania oblodzenia. Zastosowano klasyczny model matematyczny dynamiki ruchu samolotu, opisany przez układ dwunastu równań różniczkowych pierwszego rzędu. W odniesieniu do wiadomości zaczerpniętych z literatury oraz danych geometryczno- masowych samolotuTS-11 „Iskra” opracowano programkomputerowy do wyznacza- nia parametrów ruchu z uwzględnieniem oblodzenia samolotu. Na podstawie analizy wynikówprzedstawionownioski opisującewłasności dynamiczne samolotuwykonują- cego lot w warunkach oblodzenia. Manuscript received January 14, 2010; accepted for print June 15, 2011