Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 48, 2, pp. 397-413, Warsaw 2010 THE INFLUENCE OF WING TRANSVERSE VIBRATIONS ON DYNAMIC PARAMETERS OF AN AIRCRAFT Michał Wachłaczenko Air Force Institute of Technology, Aero-Engines Division, Warsaw, Poland e-mail: michal.wachlaczenko@itwl.pl The paper presents numerical analysis of airplane parameters of motion during wing transverse vibrations of an average commercial airplane. The mathematical model of the airplane motion including a system of twelve ordinary differential equations was employed. Using available pu- blications andM28”Bryza”airplanegeometrical andmass relationships, a computer code evaluating parameters of themotionwas developed. In order to take the wing transverse vibrations under consideration, the way of wing displacements in the aircraft-fixed system of coordinates was elaborated as well as vibration amplitude and frequency were as- sumed. The analysis covered two cases of wing transverse vibrations: bending of bothwing halves in phase and in counter-phase.On the basis of the obtained results, some conclusions referring the airplane dynamic properties during wing vibrations were formulated. Key words: wing vibrations, aircraft, sumulation 1. Introduction There are several forces and moments acting on the structure of an airplane during its flight in disturbed earth’s atmosphere. Because of the nature of the interactions, they can be divided into two groups. The first group of airframe loadings is connectedwith the influence of themedium surrounding and affec- ting the external surfaces of the airplane. The second one, speaking generally, are the internal interactions of the airframe and elasticity forces. They both create aeroelastic phenomena. Analysis of stability and steerability properties employs simplification that the airplane is a 6DOF (degree of freedom) rigid body.An assumption ismade that all the airframe elements occupyan invariant position toward one another throughthewhole timeof analysis.This is one sufficient condition todetermine 398 M. Wachłaczenko the mass and aerodynamic forces as well as moments acting on the aircraft during steady flight in the atmosphere. The phenomenon of aeroelasticity has accompanied the process of creation of airframes since the beginning of aviation. Varying loadings due to flow disturbances change the aerodynamic loads. This leads to vibrations of lifting surfaces (wings, tailplane) as well as bigger assemblies like the aircraft tail with the tailplane. Furthermore, the forces produced by wing vibrations are transferred through joints to the fuselage inducing its bendingor/and twisting. Aircraft vibrations cause premature wear and tear of elements and assem- blies. A proper material and structure selection allow the airframe to work below the elastic limits. Exceeding the permissible loads during flight may lead to permanent airframe deformation. Loads surpassing themaximumper- missible loading cause airframe failure and crash. Vibrations of the lifting surfaces have significant influence on the airplane dynamic parameters. For example, in the commercial aviation varying aero- dynamic forces and moments acting on the airplane, as an effect of structure vibrations, produce variable gravity loads lowering the journey comfort. Those factors worsen concentration of the air-crew and bring fatigue to passengers. The phenomenamentioned above are a kind of impulse for many research centers worldwide to make particular investigations aimed at minimising the adverse operational conditions and increasing the flight safety. Carrying out laboratory experimental research is usually insufficient, but in-flight investi- gations can be sometimes very dangerous. Therefore, parallel creation of the airframe vibration mathematical models is highly reasonable. This way, one can lower the costs of the research and the hazard of disaster. 2. Equations of motion of the airplane Spatial motion of the airplane treated as a rigid body is described by twelve nonlinear ordinary differential equations. Solving them in a general form is a complex problem which has not been accomplished so far. In practice, engi- neering calculations employ approximate methods on the basis 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 displacements of control surfaces. The influence of wing transverse vibrations... 399 The assumed simplification are: • the airplane is a rigid body with constant mass, constant moment of inertia and invariable position of the centre of mass; • control surfaces are stiff and their axes of rotation have invariable posi- tion against the airplane; • the plane Oxz is the geometrical, inertial and aerodynamic symmetry plane. The gyrostatic moment produced by elements of revolving turbo-prop en- gines is taken into consideration in the equations of motion. Figure 1 presents the systems of coordinates used in numerical analysis. Fig. 1. Systems of coordinates and angles defining their mutual position Mutual orientation of the coordinate systems are defined by the following angles: angle of attack α; angle of slideslip β; aircraft roll angle Φ; aircraft pitch angle Θ; aircraft yaw angle Ψ. Performing rotations successively by Euler angles Ψ, Θ and Φ from Oxgygzg to Oxyz, the transformation matrix is determined    x y z    =Ls/g    xg yg zg    (2.1) where 400 M. Wachłaczenko Ls/g = (2.2) =    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 from Oxayaza to Oxyz, another transformation matrix is determined    x y z    =Ls/a    xa ya za    (2.3) where Ls/a =    cosαcosβ −cosαsinβ −sinα sinβ cosβ 0 sinαcosβ −sinαsinβ cosα    (2.4) 2.1. Equations of the airplane translational motion The equation of motion for the center of mass is m dV dt =F (2.5) and can bewritten in a scalar form in the airplane-fixed system of coordinates Oxyz m(u̇+qw−rv)=X m(v̇+ru−pw)=Y m(ẇ+pv−qu)=Z (2.6) where: m is the airplanemass; V = [u,v,w]⊤ – vector of absolute velocity of the centre ofmass of the airplane; Ω= [p,q,r]⊤ – vector of rotational velocity in themobile systemof coordinates; F = [X,Y,Z]⊤ – vector of external forces acting on the airplane. Equations (2.6) are to be evaluated in the flow coordinate system Oxayaza for the sake of simplicity ofdetermination of aerodynamic forces in this system. Therefore, the velocity vector has only one component ua =V and equations (2.6) gain the following form mV̇ =Xa mraV =Ya −mqaV =Za (2.7) Assuming that the rotational velocity of the airplane-fixed system Oxyz in the inertial system is equal to Ωs and the velocity of the fixed coordinate The influence of wing transverse vibrations... 401 system Oxyz in the flow coordinate system Oxayaza is known, allows one to determine the rotational velocity of the Oxayaza system in the inertial system (Issac, 1995) Ωa =Ωs+Ωs/a =Ωs+ α̇+ β̇ (2.8) Considering that: –Ωs vector in the Oxyz system is Ωs = [p,q,r] ⊤; – β̇ vector in the Oxayaza system is β̇= [0,0,β̇] ⊤; – α̇ vector in the Oxyz system is α̇= [0,−α̇,0]⊤; andmaking use of transformations matrix (2.4), on the basis of (2.8), we get pa = pcosαcosβ+(q− α̇)sinβ+r sinαcosβ qa =−pcosαsinβ+(q− α̇)cosβ−r sinαsinβ (2.9) ra =−psinα+rcosα+ β̇ The third equation of (2.7) considers the aerodynamic force Pza as a part of the force Za depending, among others, on the rate of change of the airplane angle of attack. Thus it is calculated as follows Za =Zast+Z α̇ a α̇ (2.10) Inserting systems (2.9) and (2.10) into equations ofmotion (2.7) and trans- formaing, one obtains V̇ = 1 m (Xs cosβ+Ys sinβ) α̇= 1 cosβ− Zα̇ S mV [ZSst mV +qcosβ− (pcosα+rsinα)sinβ ] (2.11) β̇= 1 mV (Yscosβ−Xs sinβ)+psinα−rcosα 2.2. Equations of rotary motion of the airplane Thedifferential vector equation for change of the total angularmomentum is dK dt =M +Mgir (2.12) where: K = ∑ i(ri×miV i) is the total angularmomentumwith respect to the given reference point; M = [L,M,N]⊤ – vector of moment of forces acting on the airplane; Mgir = Jω×Ω – gyrostatic moment; J –moment of inertia of the engine rotor; ω= [ω,0,0]⊤ – angular velocity of the engine rotor. 402 M. Wachłaczenko Components of the gyrostatic moment, according to Kowaleczko (2003), are defined as Mgir = [Lgir,Mgir,Ngir] ⊤ = [0,−Jωr,Jωq]⊤ (2.13) As the Oxz surface is the airplane geometrical, inertial and aerodynamic symmetry plane and the engine rotors revolve with the angular velocity ω along an axis parallel to the Ox axis, the scalar form of equation (2.12) 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 (2.14) Izṙ− (Ix− Iy)pq− Ixz(ṗ− qr)=N+Jωq Transforming system (2.14) ṗ= 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)] (2.15) ṙ= 1 IXIZ − I 2 XZ {[L+(IY − IZ)qr+ IXZpq]IXZ + +[N+Jωq+(IX − IY )pq− IXZqr]IX} 2.3. Modified equations of motion used in numerical analysis To equations (2.11) and (2.15) we need to add equations allowing for de- finition of the current Euler angles (kinematic relations) of the airplane Φ̇= p+(rcosΦ+q sinΦ)tanΘ Θ̇= qcosΦ−r sinΦ (2.16) Ψ̇ =(rcosΦ+q sinΦ) 1 cosΘ and equations determining the position of the airplane centre of mass in the gravitational system of coordinates Oxgygzg ẋ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Ψ)] (2.17) żg =V [−cosαcosβ sinΘ+sinβ sinΦcosΘ+sinαcosβ cosΦcosΘ] The influence of wing transverse vibrations... 403 Equations (2.11) and (2.15)-(2.17) create a set of twelve nonlinear ordina- ry differential equations describing spatial motion of the airplane treated as a rigid body with constant mass. The vector of motion parameters has the following components: V , α, β, p, q, r,Φ,Θ, Ψ, xg, yg, zg. 3. Forces and moments acting on the airplane The resultant vector of external forces which affect the airplane during flight (right-hand side of equation (2.5)) is a sum of the following forces F =Q+T +R (3.1) where: T is the thrust force; R – aerodynamic force; Q – terrestrial gravity force. The corresponding components of resultant external force (3.1) can be written as Xa =Qxa+Txa+Rxa Ya =Qya+Tya+Rya Za =Qza+Tza+Rza (3.2) 3.1. Gravity force Q The terrestrial gravity force Q acting on the airplane in the Oxgygzg system has only one component Q = [0,0,mg]⊤. Using relations (2.2) and (2.4), these components in the Oxayaza system can be determined    Qxa Qya Qza    =L−1 s/a Ls/g    0 0 mg    = (3.3) =    mg(−cosαcosβ sinΘ+sinβ sinΦcosΘ+sinαcosβcosΦcosΘ) mg(cosαsinβ sinΘ+cosβ sinΦcosΘ− sinαsinβcosΦcosΘ) mg(sinαsinΘ+cosαcosΦcosΘ)    3.2. Engine thrust force T Both engine propeller thrust vectors TL and TR lie in planes parallel and equidistant to the aircraft symmetry plane Oxz as well as parallel to the Ox 404 M. Wachłaczenko axis. In the Oxyz system, T has only one component T = [TL +TP ,0,0] ⊤. Using transformation matrix (2.4), the thrust force vector is obtained    Txa Tya Tza    =L−1 s/a    TL+TP 0 0    =    (TL+TP)cosαcosβ −(TL+TP)cosαsinβ −(TL+TP)sinα    (3.4) 3.3. Aerodynamic force R Projections of the resultant aerodynamic R force on axes of the Oxayaza system are Rxa =−Pxa =−Cxa ρV 2 2 S Rya =−Pya =−Cya ρV 2 2 S (3.5) Rza =−Pza =−Cza ρV 2 2 S where: Cxa,Cya,Cza denote the drag, side, lift coefficients, respectively; ρ is the air density; S – wing (lifting surface) area. Fig. 2. Forces acting on the airplane Themoment M = [L,M,N]⊤ on the right-hand side of equation (2.12) is the resultantmoment affecting the airplane, whose components are defined as L=Cl ρV 2 2 Sl M =Cm ρV 2 2 SbA N =Cn ρV 2 2 Sl (3.6) The influence of wing transverse vibrations... 405 where: Cl, Cm, Cn are the roll, pitch, yaw moment coefficients, respectively; l – wing span; bA – mean wing chord. 3.4. Coefficients of the aerodynamic force and moment In this paper, experimentally obtained static aerodynamic characteristics Cxa(α,Ma), Cza(α,Ma),Cm(α,Ma) are used. Along with the determined de- rivatives they allow one to define final expressions defining force andmoment coefficients: – drag force: Cxa =Cxa(α,Ma) – side force: Cya =C β yaβ+C p yap+C r yar+C δV ya δV +C δl yaδl – lift force: Cza =Cza(α,Ma)+C α̇ zaα̇+C q zaq+C δH za δH – rolling moment: Cl =C β l β+C p l p+C r l 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 3.5. Initial conditions of flight Originally, the airplane performs a rectilinear steady flight with the given constant speed. For this condition, we can take: p= q= r=0 β=0 Φ=0 γ =Θ−α=0 V̇ = α̇= β̇= ṗ= q̇= ṙ= Θ̇= Φ̇=0 In order to satisfy these conditions, the resultant foce andmoment acting on the airplanemust be equal to zero Pxa sinα+(Pza−mg)cosα=0 (3.7) T −Pxacosα+(Pza−mg)sinα=0 After necessary transformations, equations for the steady flight angle of attack α (3.7)1 and indispensable thrust T (3.7)2 are obtained. Knowing the airplane angle of attack α, we determine the elevator displacement angle δH0, so the coefficient of pitching moment is equal to zero. 4. Simulation of wing transverse vibrations In simulation ot thewing vibration, themodel ofM28 ”Bryza”Polish airplane with strut-bracedwingswas used.Along the span, thewing has varying cross- section areas (different moments of inertia), mass (different mass loads) as 406 M. Wachłaczenko well as varying aerodynamic loads. Some assumptions were made in order to simplify the model of wing vibrations such as: • the wing has a constant cross-section area along its span; • the wing is a fixed-free beam supported at yL (yP) from the wing root (airplane symmetry plane); • thewing operates under a continuous aerodynamic load, invariant along the span; • moments of inertia and moments of deviation are constant through the simulation. On the basis of the above assumptions, the bending moment causes its parabolic-like deflection (Fig.3.). Because of the scope of this paper, the pro- blem of finding the frequency of structural vibrations and their amplitude has not been developed. Fig. 3. Projection of the deflected wing (along the flight direction) In the process of numerical simulation, the wing is divided into 50 evenly distributed elements. Because of thewing taper and sweep, each element has a different area and point of application of the aerodynamic force. The aerody- namic force is applied at the 25% chord of the considered wing element. The distance between the centre of the wing element and the airplane symmetry planemay be obtained according the formula below ylok(ne)=± (∆l 2 +(ne−1)∆l ) (4.1) where ”+” refers to the right half of the wing, and the ”−” to the left half; ∆l is thewidth (span) of thewing element; ne – number of thewing element. Considering a given element of the vibrating wing, one needs to determine the local angle of attack (further named: AOA). Assuming that the airfoil chord is parallel to the airplane longitudinal axis Ox, the total angle of attack The influence of wing transverse vibrations... 407 of the wing element (according to Fig.4.) is equal to the sum of the airplane AOA and the change of AOA generated by thewingverticalmovement during transverse vibrations. The above may be written as αi(t)=α0(t)±∆αi(t) ±∆αi(t)= arcsin (±Vzi(t) V (t) ) (4.2) where: ”+” indicates that the elementmoves ”downward”, and”−””upward”; α0 is the airplane angle of attack; ∆αi – change of AOA produced by trans- verse vibrations; αi – current AOA value of the element. Fig. 4. Vertical velocity and AOA of the wing element Fig. 5.Wing tip displacement and vertical velocity during transverse vibrations The component calculated in equation (4.2)2 produces the damping force for theverticalmovement,which leads toadecrease in amplitudeof vibrations. The lift force for eachwing element is derived.Calculations of the lift force take into consideration the wing sweep, local AOA value and distribution of circulation over the full span of the lifting surface. Figure 6 shows decomposi- tion of the lift force along the wing span for aircraft velocity of 70m/s. 408 M. Wachłaczenko Fig. 6. Distribution of the wing lift force: 1 – elementary lift forces not including circulation; 2 – elementary lift forces including circulation; 3 – elementary lift forces for wingmoving ”downwards”; 4 – elementary lift forces for wingmoving ”upwards” The total lift force produced by the wing of the aircraft is equal to P totalza = l/2 ∫ −l/2 Pza = 0 ∫ −l/2 Pza+ l/2 ∫ 0 Pza =P L za+P P za For determination of the lift force used in the consideredmodel, the men- tioned force may be obtained from equation (4.3) P totalza =P L za+P P za = 25 ∑ i=1 (PLza)i+ 25 ∑ i=1 (PPza)i (4.3) The wing oscillates with a given frequency and amplitude. In each ite- ration of the numerical simulation, the vertical component zel of the wing element position (Fig.3.) and vertical velocity Vz are evaluated. Figure 5 pre- sents trajectory of the wing tip and its normal velocity. For each element, the coefficients of drag, lift and pitching moment are determined. 5. Numerical analysis of perturbed airplane motion In order to analyse the influence of wing transverse vibrations, using a com- puter code for solving the equations ofmotion,many cases of vibrations of the lift surfaces may be simulated. Because of the complex structure of the wing causing difficulties in estimating the vibration frequencies and amplitudes and the lack of airplane in-flight experimental data, the mentioned above values were evaluated (Issac, 1995). The influence of wing transverse vibrations... 409 Both cases of airplanemotionwithwingvibrations are basedon themathe- maticalmodel described in Section 4. The steering vector XST = [δH,δV ,δl] ⊤ remains invariant in the analysis, which corresponds to a non-controlled flight condition. The airplane performs a steady flight in smooth configuration (flaps on δKL = 0 ◦) with velocity of V = 70m/s (∼ 252km/h) at altitude of H =1000m. The initial turboprop engines thrust in the analysed time interval rema- ins constant and has been determined on the base of steady flight condition (Eq. (3.7)2). The static pitching moment is compensated by deflection of the horizontal control surface, which depends on the flight speed and is equal δH0 =−0.44 ◦. The system of twelve nonlinear ordinary differential equations (2.11) and (2.15)-(2.17) was used to carry out the analysis. Figures 7-14 present the results of airplane motion simulation with wing transverse vibrations covering two cases: (a) – ”in phase” wing vibrations – bothwing halves aremoving in the same direction (along the airplane normal axis Oz); (b) – ”in counter-phase” wing vibrations – wing halves are moving in the opposite direction. Calculations were carried out for 300 seconds. Thewing vibration occured after 10th seconds of simulation and disappeared at the 20th second. There was no pilot response, so the behaviour of the airplane itself was investigated. Fig. 7. Airplane velocity Concerning thecase (a), there is a slight loss in theairplanevelocity (Fig.7) and an increase in the angle of attack (Fig.8) as the wing tips move upward at the very first moments of simulation. There is a temporary increase in the pitching moment resulting in the nose ”pull up” (Fig.10). Because both wing halves move in the same direction, the side forces produced by the lift force are equal and opposite. Thus a small side-slip angle is induced by the gyrostaticmoment fromengines rotors.Because of this, the roll andyawangles 410 M. Wachłaczenko Fig. 8. Airplane angle of attack Fig. 9. Roll angle Fig. 10. Pitch angle Fig. 11. Yaw angle The influence of wing transverse vibrations... 411 Fig. 12. Flight altitude Fig. 13. Inclination angle Fig. 14. Airplane trajectory (Fig.9 andFig.11) are small too (coupling of longitudinal and lateralmotion). A decrease in the altitude (Fig.12) and a small deviation from the initial direction of flight is observed (Fig.14). In case (b),when thewinghalvesmove in the opposite direction, the above mentioned elementary side forces produce bigger side-slip angles as well as the roll and yaw. This leads to increasing deviation from the initial direction of flight. In comparison, there is amajor velocity increase and a decrrease in the angle of attack at the very first moments of simulation than it was observed in case (a). 412 M. Wachłaczenko 6. Conclusions On the basis of results of numerical simulation, it is possible to state that a disturbed numerical model of M28 ”Bryza” airplane remains stable because the parameters of motion return to their primary values: • longitudinal stability → velocity, AOA and pitch angle return to the fixed values; • lateral stability → roll angle returns to 0 and yaw angle is stabilized at a certain level. Comparing two cases of wing transverse vibrations, it can be found that the ”in counter-phase” wing vibrations have considerable influence on flight parameters. Because of the coupling of longitudinal and lateral motions, the growing side-slip angle produces bigger variations of the parameters. In both cases, a little drop of the altitude is noticed. During wing vibrations, all mo- tion parameters oscillate, which results in unfavourable working conditions for many onboard systems, for instance the autopilot system. The vibrations badly influence the working environment for crew and travel comfort as well. Wing vibrations have also significant effect on the airframe loading and durability. Varying loads lead to premature tear andwear of wing joints in the fuselage, suspensed devices andmountings of the control surfaces in the wing. Deflections of the wing along the airplane normal axis Oz can produce additional forces acting onailerons,whichmay result inbiggerwingdeflections and loads in the lateral control system. References 1. Fiszdon W., 1962,Mechanika lotu, Cz. I, II, PWN,Warszawa 2. Hodges D.H., Pierce G.A., 2002, Introduction to Structural Dynamics and Aeroelasticity, Cambridge Aerospace Series, Cambridge 3. Issac J.C., 1995, Sensitivity Analysis of Wing Aeroelastic Responses, Virginia Polytechnic Institute, Blacksburg 4. Kowaleczko G., 2003, Zagadnienie odwrotne w dynamice lotu statków po- wietrznych, WydawnictwoWAT,Warszawa 5. Krzyżanowski A., 1984,Mechanika lotu, skryptWAT,Warszawa 6. Ostosławskij N.W., 1957,Aerodinamika samoleta, Oborongizdat,Moskwa The influence of wing transverse vibrations... 413 7. Szabelski K., Jancelewicz B., 2002, Wstęp do konstrukcji śmigłowców, WKŁ,Warszawa 8. QinZ., 2001,Vibration andAeroelasticity of AdvancedAircraftWingsModeled as Thin-Walled Beams, Virginia Polytechnic Institute, Blacksburg Wpływ drgań giętnych skrzydła na dynamiczne parametry lotu samolotu Streszczenie W pracy przedstawiono wyniki numerycznej symulacji drgań giętnych skrzydła samolotu dyspozycyjnego. Zastosowanoklasycznymodelmatematyczny dynamiki ru- chu samolotu, opisanyprzezukładdwunastu równańróżniczkowychpierwszego rzędu. W celu uwzględnienia drgań skrzydła zamodelowano sposób przemieszczania się ele- mentówskrzydław samolotowymukładziewspółrzędnychoraz założonoczęstotliwość i amplitudę drgań skrzydła.Analizę przeprowadzono dla dwóch przypadków zginania skrzydła: zginanie połówek skrzydła w fazie i w przeciwfazie. Na podstawie analizy opracowanownioski dotyczącewłasności dynamicznych samolotu podczas drgań jego powierzchni nośnych. Manuscript received September 23, 2009; accepted for print October 14, 2009