Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 41, 4, pp. 823-852, Warsaw 2003 A STUDY OF THE FLIGHT DYNAMICS HELICOPTER CARRYING AN EXTERNAL LOAD USING BIFURCATION THEORY AND CONTINUATION METHODS Krzysztof Sibilski Department of Aeronautical Engineering, Wrocław University of Technology and Air Force Institute of Technology, Warsaw e-mail: krzysztof.sibilski@pwr.wroc.pl The paper presents a study of the flight dynamics of a helicopter with an articulated rotor, carrying a suspended load. The aircraft model includes rigid body dynamics, individual flap and lag blade dynamics, and inflow dynamics. The load is a point mass with a single suspension point. Results are obtained for load masses of up to 1500kg, with load-to-helicopter mass ratios up to 33%, and cable lengths up to 55m. The presence of the external loadmodifies the flight dynamics and handling quality characteristics of the helicopter because the dynamic and aerodynamic characteristics of the load maymake it unstable in certainfight conditions.Maneuvers of this system in wide flight regimes involves non-linear aerodynamics and inertial coupling. In can be stated that the helicopter with the suspended load is a inherently non-linear and time varying system. Theory of dynamical systems provides amethodology for studying such non-linear systems. Bifurcation theory is a part of that theory. It considers changes in the stability of the systemwhich lead to qualitatively different responses of it. In this paper, results from the theory of dynamical systems are used to predict the nature of instabilities causedbybifurcations and the response of the rotorcraftwith the suspended load that follow such bifurcations. Key words: non-linear helicopter dynamics, bifurcation theory, continuation method 1. Introduction Application of helicopters to transport of heavy and bulky loads creates stability problems, especially in the hovering with hanging loads. Bothmilita- ry and commercial operators have exploited the capability of a helicopter to 824 K.Sibilski rapidly move heavy loads to locations where the use of ground based equip- ment would be impractical or impossible. The external load can modify the flight dynamic characteristics of a helicopter because the load behaves like a pendulum,and can change the natural frequencies andmode shapes of low fre- quencymodes of the helicopter. Also, the aerodynamics of the loadmaymake it unstable in certain flight conditions, with repercussions on the stability and safety of the entire helicopter/load system. The dynamics of a helicopter with external suspended loads received con- siderable attention in the late 1960’s and early 1970’s. Therewere two reasons for that interest. Firstly, the extensive external load operations in the Viet- namwar, and secondly, the Heavy-Lift Helicopter program (HLH). This inte- rest has been renewed recently, prompted by the re-evaluation and extension of the ADS-33 [1] Helicopter Handling Qualities Specifications to transport helicopters, and by the expectation of emerging new cargo helicopter manu- facturers. 1.1. Historical outlook One of the first theoretical studies of the dynamics of a helicopter with a slung load is due to Lucassen and Sterk (1965). A simple 3-degree of freedom model of the hover longitudinal dynamics of the helicopter and the angular displacement of the loadwas used.A single suspensionpointwas assumed and the aerodynamic forces and moments on the load were neglected. In general, the pole associated with the load pendulummode was stable; the phugoid re- mained unstable, but its frequency decreasedwith the increasing cable length. For some combinations of parameters, the helicopter mode became unstable while the load mode was stabilised. Szustak and Jenney (1971) pointed out that a conventional stability augmentation system was not adequate for pre- cising the hover and load release, and could result in pilot-induced oscillations (PIO). A more effective solution consisted of an inner loop, in which the re- lative motion of the aircraft and load was fed back to cyclic, and an outer loop, in which the aircraft position above the ground was fed back, again to cyclic. Dukes (1973a) studied the basic stability characteristics of a helicopter with a slung load, andpossible feedback stabilisation schemes, andappropriate piloting strategies for various manoeuvres (Dukes, 1973b). In this work a 3-degree of freedom longitudinal helicopter/load model was used. The positive pitch damping, whether provided by the rotor alone or also by the flight control system, did not necessarily increase the stability of the pendulum mode of the load. This mode, essentially undamped, could becomeunstable for certain configurations (i.e. cable lengths, loadweight, and A study of the flight dynamics helicopter... 825 relative position of the attachment point and gravity center of the aircraft). Thepitch dampingprovided at best amoderate increase in the dampingof the mode. A feedback control scheme in which the attachment point was actively and longitudinally moved, proved very effective in the paper, but its practical feasibility was not explored. The previous studies were limited to hover or low speed flight, and therefore the aerodynamics of the suspended load did not play a significant role. Slung loads are rarely aerodynamically shaped bodies. Typical loads are bluffbodies thatmaybe subject todynamic instabilities triggered byunsteady aerodynamics. Poli andCromack (1973) studied the stability in forward flight of a helicopter carrying a container and a circular cylinder. The results indica- ted that long cables, high speeds, and lowweights increased the stability of the loads. A stability study in forward flight by Cliff and Bailey (1975) partially confirmed the results byPioli andCromack (1973) because the decrease of the weight improved the stability, but longer cables were found to be destabili- sing. The differences may be due to different aerodynamics of the load, which wasmuchmore idealised by Cliff and Bailey (1975). The lowering of the drag increased the stability. The lateral and longitudinal stability was governed by the same parameters, but the conditions for the lateral stability proved more stringent. The results presented byŁucjanek andSibilski (1978) confirmed the result obtained by Cliff and Bailey. A few years later, Sibilski and Łucjanek (1983) addressed the stability of a single-point configuration of the load su- spension system. The analysed model was much more sophisticated than in anyof thepreviouslymentioned studies, and included full non-linear equations for the helicopter motion and dynamics of the 3-degrees of freedom model of the suspended load. The equations were then linearized for the stability ana- lysis, and the effect of several configuration parameters was investigated. The cable length, fore/aft and vertical position of the suspension point, and load weight were all found to affect the stability. Depending on the combination of parameters, somemodes could be stabilised and others destabilised, but in all, the instabilitieswerequiteweak.Concurrently,Nagabhushan (1985) addressed in his analysis the low-speed stability of a single-point configurationof the load suspension system. The analysed model included full non-linear equations for the rigid body aircraft motion and rotor flap dynamics. His results confirmed that the cable length, fore/aft and vertical position of the suspension point as well as the load weight affect the stability of the helicopter-suspended load system. More recently, Cicolani et al. (1989) reported the results of flight tests of a UH-60 helicopter, including frequency responses obtained using system 826 K.Sibilski identification techniques. While the study focused primarily on the system identification and simulation validation, several conclusions were presented regarding the effect of the loads on the flight dynamics and handling quality characteristics. The increasing load weight reduced the lateral bandwidth; a further increase could reduce the bandwidth to a value below that of the pendulumfrequency.The longitudinal stabilitymarginswerenotvery sensitive to the load, but the lateral stability margins were degraded. The effect on the bandwidth and phase delay was highly variable depending on the load configuration.At last, Fusato et al. (1999) explored some fundamental aspects of the dynamics of a helicopter with an articulated rotor with an external load suspended from a single attachment point. The results indicated that the external load affects the trimstate primarily through theoverall increase in the weight of the aircraft, both in straight and in turning flight. The influence of the length of the cable was negligible. Substantial dynamic coupling with the Dutch roll mode occurred. The load mode consisted primarily of the lateral motion of the load. The effect of the load on the phugoid was very small. The suspended loadmodified the on-axis roll frequency response by adding a notch to the gain curves and a roughly 180-degree jump in the phase curves. The modifications of the frequency response introduced by the load occurred primarily at frequencies lower than thoseused todetermine thebandwidthand phase delay according to the ADS-33 specifications. The phase shifts caused additional crossings of the 135-degree delay line, which at least formally, can reduce the phase bandwidth considerably. If these additional crossings were ignored, the changes in the bandwidth and phase delay would be generally small. 1.2. Subject of the present investigations Ahelicopter carrying a suspended load is a inherently non-linear and time varying system. Therefore, a linear model is adequate for basic studies of the flight dynamics of helicopters with suspended loads, but it cannot describe some important practical problems. For example, it is not possible to model, using a linear approach, the sling load ”vertical bounce” phenomenon. Ano- ther type of problems that linear models cannot capture is the aerodynamic instability due to unsteady inflow, and/or the non-streamlined shape of many suspended loads. These instabilities, described by Gabel and Wilson (1968), Poli and Cromack (1973), Sheldon (1977), Simpson and Flower (1981), and Cicolani et al. (1995) often limit the maximum speed of helicopters. Recent developments in the field of numerical analysis of non-linear equ- ations created a class of computer algorithms known as the continuation me- A study of the flight dynamics helicopter... 827 thods. Those methods use predictor- corrector techniques to follow solution curves of systems of non-linear equations of motion represented by functions of a number of variables and parameters, respectively. This approachwas suc- cessfully demonstrated for the flight dynamics analysis a helicopter with a suspended load. The continuation methods are a class of predictor corrector techniques for the solution of systems of non-linear algebraic equations, which are functions of a number of parameters, over a specified range of the pa- rameters. The general technique is to fix all parameters but one, and trace the steady states of the system as a function of this parameter. The stabili- ty of each steady state can be determined by calculating the eigenvalues of the linearized system. Any changes in the stability from one steady state to the next will signify bifurcation. Theory dynamical systems has provided a powerful tool for analysis of non-linear phenomena of the aircraft behaviour. In the application of this theory, the numerical continuation methods and bi- furcation theory have been used to study roll-coupling instabilities, stall/spin phenomena, and dynamics of high angles of attack is a number of aircraftmo- dels. Results of great interest have been reported in several papers (one should mention here papers byCarroll andMehra (1982), Guicheteau (1990), Jahnke and Culick (1994), Avanzini and de Matteis (1998), Sibilski (1999a,b, 2000), Marusak et al. (2000)). The continuation methods are numerical techniques for calculating the steady states of systems of ordinary differential equations, and canbeused to study the roll coupling instabilities andhigh-angle of attack instabilities. The objectives of this paper are: • To present formulations and solutions of a mathematical model of a helicopter with an articulated rotor, carrying an external load • To study the effects of the cable length and load weight on the helicop- ter/load system dynamic characteristics. In the present study some fundamental aspects of the nonlinear dynamics of a helicopter with an articulated rotor and with an external suspended lo- ad are studied using the continuation and bifurcation methods, by means of checking the stability characteristics related to unstable equilibria. Numerical simulations are used to verify the predictions. The mathematical model of a helicopter used in this study is a nonlinear blade element type model that includes the fuselage, rotor, main rotor inflow, and propulsion system dyna- mics. The 6-degrees of freedom rigid bodymotion of the aircraft, 2-degrees of freedom rigid bodymotion of the slung load, and an articulated, four-bladed main rotorwith rigid blades are assumed.Theaerodynamic loads areunsteady forces and moments in the direction determined by the local airflow (defined 828 K.Sibilski by the suspended load angle of attack and slip angle). Unsteady aerodynamic effects aremodeled using theONERAdynamic inflowmodel (Tran andPetot, 1981). The state vector has a total of 30 elements: flap and lag displacements and rates for each of the 4 blades (16 states); 9 rigid body velocities, rates, and attitudes; and 2 external load angles with their respective rates. The for- malism of analytical mechanics allows one to present dynamic equations of motion of the coupled load/helicopter dynamic system giving an incredibly interesting and comfortable tool for the construction of equations of motion. An example can be the Boltzmann-Hamel equations, which are a generaliza- tion of the Lagrange equations of the second kind in quasi-coordinates. These equations are written in the form allowing one to create proceduresmeant for their automatic formulation, (e.g., by means of the well known commercial software asMathematica R© orMaple V R©). 2. Mathematical model of the helicopter – slung load system 2.1. Systems of coordinates The external load is essentially modeled as a point mass that behaves like a spherical pendulum suspended from a single point. The geometry and the relevant coordinate systems are shown in Fig.1. The position of the load is described by the three angles ΘL, ΨL, and ΦL, where ΨL is the azimuth angle of the load, ΦL is measured from the y ′ 3 axis, and ΘL is measured from the z4 axis. The following systems of coordinates (Fig.1) are used in futher constri- butions: – Averticalmoving systemof co-ordinates Oxgygzg, the Ozg axis ofwhich is vertical and directed downwards. – A system of co-ordinates Oxyz attached to the rotorcraft. The origin of this system overlap the rotorcraft centre of mass. The Oxz plane coincides with the symmetry plane of the aircraft, and the Oz axis is directed downwards. – A system of co-ordinates O4x4y4z4, with the origin O4 that overlaps the suspension point. All axes are parallel to the slung load axis of inertia, the Ox4 axis of that system is directed forwards, and the Oz4 axis is directed downwards. – A system of co-ordinates O5x5y5z5 attached to the suspended load. The origin of this system overlap the suspended load center of mass. The A study of the flight dynamics helicopter... 829 Fig. 1. Systems of co-ordinates axes of this system are parallel to the axis of the system of co-ordinates Oxyz. – A system of co-ordinates O6x6y6z6 attached to the suspended load. The origin of this system overlap the suspended load center of mass. The axes of this system are parallel to the axis of the system of co-ordinates O4x4y4z4 (the axes of this system overlap the slung load axis of inertia). The systems of coordinates attached to the main rotor blades are shown in Fig.2. The following systems of coordinates are used: – A system of co-ordinates OxWNyWNzWN with the origin that overlap the centre of the main rotor hub. All axes are parallel to the system of co-ordinates Oxyz attached to the aircraft. – Systems of co-ordinates of themain rotor hub OrisizWN (i =1,2, ...,n; n – number ofmain rotor blades) attached to themain rotor hub.Those systems turn with the rate Ω of the main. 830 K.Sibilski Fig. 2. Systems of co-ordinates attached to the main rotor andmain rotor blades – Systems of co-ordinates attached to the main rotor blades. The rotor blade is mounted to the hub on a universal join – free to flab (flapping hinge PW , the systemof co-ordinates PWtβbβnβ), leador lag (laghinge PO, the systemof co-ordinates POtζbζnζ), butfixed inpitch (feathering hinge, the system of co-ordinates POtbn). 2.2. Equations of motion The formalism of theoretical mechanics allows one to present dynamic equations of motion of mechanical systems in quasi-coordinates, giving an in- credibly interesting and comfortable tool for the construction of equations of motion of an aircraft. An example can be the Boltzmann-Hamel equations, whichare a generalisation of theLagrange equations of the secondkind expres- sed in quasi-coordinates. The Boltzmann-Hamel equations have the following form (Gutowski, 1972; Arczewski and Pietrucha, 1993) d dt (∂T∗ ∂ωσ ) − ∂T∗ ∂πσ + k ∑ µ=1 k ∑ λ=1 γ µ σλ ∂T∗ ∂ωµ ωλ = Q ∗ σ (2.1) where T∗ is the kinetic energy (function of quasi-coordinates), ωσ – quasi-velocity, πσ – quasi co-ordinate, qλ,qσ – generalized co-ordinates, Q∗σ = ∑k σ=1Qσbσµ – component of the generalized force vector, k – num- A study of the flight dynamics helicopter... 831 ber of degrees of freedomof themechanical system, γrµα –Boltzmann symbols (Gutowski, 1972; Arczewski and Pietrucha, 1993) γrµα = k ∑ r=1 k ∑ α=1 (∂arσ ∂qλ − ∂arλ ∂qσ ) bσµbλα (2.2) where arσ, brσ are elements of the transformation matrix. Relations between the quasi-velocities and generalized velocities are shown in the equations ωσ = k ∑ α=1 aσα(q1,q2, ...,qk)q̇α (2.3) q̇σ = k ∑ m=1 bσµ(q1,q2, ...,qk)ωµ σ =1, ...,k Equations (2.3) can be written in the matrix form Ω=AT q̇ q̇=A −1 T Ω=BTΩ (2.4) where Ω is the vector of the quasi-velocities and q – vector of the generalized co-ordinates Ω= [ω1,ω2, ...,ωk] ⊤ (2.5) q= [q1,q2, ...,qk] ⊤ The construction of thematrix AT depends on the explored issue. In the case whenwe consider amodel of a helicopter carrying a suspended load treated as a system containing a rigid fuselage and n rigid blades of themain rotor fixed to thehubbymeans of threehinges, and threedegrees of freedomhinging load, the quasi-velocities and generalized co-ordinates have the following forms Ω= [u,v,w,p,q,r,Ω,β̇1, ..., β̇n, ζ̇1, ..., ζ̇n, Ψ̇L,Θ̇L, Φ̇L] ⊤ (2.6) q= [xs,ys,zs,Φ,Θ,Ψ,ψ,β1, ...,βn,ζ1, ...,ζn,ΨL,ΘL,ΦL] ⊤ Thematrix AT is then as follows AT =    AG 0 0 0 CT 0 0 0 I    (2.7) 832 K.Sibilski Thematrices AG and CT are classicalmatrices of transformationsofkinematic relations, and can be foundbySibilski (1998), the unitmatrix Ihas dimension (3n+1)× (3n+1), n is the number of the main rotor blades. Thematrices Di can be determined as follows Di = dai dq =        ∂a11 ∂q1 · · · ∂a1k ∂qk · · · · · · · · · ∂ak1 ∂q1 · · · ∂akk ∂qk        (2.8) where the vector ai means ith row of the matrix AT . In the matrix notation, the Boltzmann symbols can be presented in the form of elements of the block matrix Γ(k×(k×k)) with k being the number of degrees of freedom of the examined dynamic system Γ=      Γ1 Γ2 ... Γk      =      B ⊤ T (D1−D ⊤ 1 )BT B ⊤ T (D2−D ⊤ 2 )BT ... B ⊤ T (Dk−D ⊤ k )BT      (2.9) Thematrix Γ can be presented in the short form Γ=B⊤T (D−D ⊤)BT (2.10) Finally, the Boltzmann-Hamel equations written in the matrix form can be presented as follows d dt (∂T∗ ∂Ω ) +(Γ⊤Ω) ∂T∗ ∂Ω −B ⊤ T ∂T∗ ∂q =Q−U (2.11) The vector Q is the sum of aerodynamic loads and another non-potential forces acting on the helicopter-slung load system, U is the vector of potential forces acting on the helicopter and external hanging load. Equation (2.11) is very comfortable to use in procedures of automatic formulation of equations of motion. In our case, the subject of considerations is the problem of dynamics of a helicopter with a slung heavy load. The quasi-velocities vector is given by Eq. (2.6)1. The total kinetic energy of the system is the sum of the kinetic energy of the rigid fuselage of the helicopter, rotor blades, and the slung load T∗ = T∗H + n ∑ i=1 T∗Bi+T ∗ L (2.12) A study of the flight dynamics helicopter... 833 According to the general theorem, the kinetic energy of the airframe is (see Sibilski, 1998) T∗H = 1 2 mV 2+ 1 2 Ω ⊤ kJAΩk (2.13) The kinetic energy of the ith rotor blade is (see Sibilski, 1998) T∗Bi = 1 2 mBi[V+Ωk×xBi+(Ωk+ωBi)×RBi] 2+ 1 2 (Ωk+ωBi) ⊤ JBi(Ωk+ωBi) (2.14) where xBi is the vector connecting the centre of gravity of the aircraft with the centre of the main rotor hub, RBi – vector connecting the centre of the main rotor hubwith the blade centre of gravity, JBi –moment of inertia of the ith rotor blade with respect to the flapping axis, ωBi – vector of the relative angular velocity of the rotor blade ωBi = Ωzwn+ dβi dt tβi + dζi dt nζi + dθi dt b (2.15) The kinetic energy of the slung load is given by the following formula T∗L = 1 2 mL[V+Ωk×H+(Ωk+ωL)×RC] 2+ 1 2 (Ωk+ωL) ⊤ IL(Ωk+ωL) (2.16) where H is the vector connecting the helicopter centre of gravity with the point of suspension of the external load, RC – vector connecting the point of suspension with the centre of gravity of the external slung load, IL – tensor of inertia of the slung load. Aftermaking some transformations, the relation for the kinetic energy can be presented in the form T∗ = 1 2 Ω ⊤ EΩ (2.17) Thematrix Edepends on themass distribution of the airframeand control surfaces, and has the form E=     M −S1 S (E) 2 S JA N (E) (S (E) 2 ) ⊤ (N(E))⊤ I (E) S     (2.18) Elements of the matrices M, N(E), JA, S1, S (E) 2 , I (E) L , can be found in Sibilski (1980, 1998). 834 K.Sibilski After diferentiation and transformations, we obtain a set of equations de- scribingmotionof thehelicopterwithanarticulated rotor, carryinganexternal suspended load EΩ̇ + [ (Γ⊤Ω)E−B⊤Ω⊤ dE dq ] Ω=Q−U (2.19) Equation (2.19), together with the kinematic relations, create a set of non- linear ordinary differential equations of the first order that describemotion of the helicopter-suspended load system. Thekinematic relations canbe found inEtkin (1972), Sibilski andŁucjanek (1983), Sibilski (1980, 1998). These equations arewritten in the form allowing onetocreateproceduresmeant for their automatic formulation, (e.g., bymeans of the well known commercial software as Mathematica R© or Mathcad R©). 2.3. Modeling of aerodynamic loads Precise description of aerodynamic forces and moments found in the equ- ations of motion is the fundamental source of difficulties. In each phase of the flight, the dynamics and aerodynamics influence each other, which disturbs the precise mathematical description of these processes. The requirement for a method enabling calculation of the aerodynamic load stem both from the flowenvironmentand fromthealgorithmsused in theanalysis of thehelicopter flight. Thebifurcation approach is very fruitfulwhen the sources andnature of aerodynamic phenomena are considered. It is assumed that the airframe mo- del consists of the fuselage, horizontal tail, vertical tail, and landing gear. The fuselagemodel is based onwind tunnel test data (as a function of the angle of attack α and the slip angle β). The horizontal and vertical tails are treated as aerodynamic lifting surfaces with lift and drag coefficients computed from data tables as functions of the angle of attack α and the slip angle β. The tail rotor is a linearmodel using the strip-momentum theorywith a uniformly distributed inflow. The influence of the rotor dynamic inflow on the airframe and aerodynamic forces and moments of the suspended load are included in the model. The technique used provides the essential effects of the increased interference velocity with an increased rotor load and decreased interference as the rotorwake deflects backwardwith the increased forward speed, seeNar- kiewicz (1994). On the basis of the results presented byPit andPeters (1981), Krotophalli et al. (1999), Ermentrout (2001), the effects of changing velocity due to helicopter angular rates are included. Special techniques are proposed to calculate the aerodynamic forces andmoments acting on the external load taking into account the dynamics inflow and the interference of the rotor wake A study of the flight dynamics helicopter... 835 on the slung load local inflow velocity and angles of attack and slip (Sibilski, 1980). The aerodynamic data for a NACA 23012 airfoil in the range of the angle of attack ±23◦, and the compressibility effects were included. Those data were blended with suitable low speed data for the remainder of the 360◦ range to model the reversed flow region and fully stalled retreating blades. Semi-empirical methods, that use differential equations, were used to predict the unsteady aerodynamic loads and dynamic stall effects. The model deve- loped by ONERA (Tran and Petot, 1981) was used. The ONERAmodel is a semi-empirical, unsteady, non-linear model, which uses experimental data to predict aerodynamic forces on an oscillating airfoil which experiences dynamic stall. State variable formulations of aerodynamic loads allow one to use the existing codes for rotorcraft flight simulation. 3. Dynamical systems theory In this paper, we will study equations of the following form ẋ = f(x,t;µ) (3.1) and x 7→ g(x;µ) (3.2) with x ∈ U ⊂ ℜn, t ∈ ℜ1, and µ ∈ V ⊂ ℜp, where U and V are open sets in ℜn and ℜp, respectively. We view the variables x as a vector of n state variables, the variables µ as a vector of m parameters (or controls), ẋ is the time derivative of x, and f : ℜn ×ℜm → ℜn is the smooth vector field (n non-linear functions). Note that both open loop (uncontrolled) and closed loop rigid-body flight dynamical systems can usually be represented in the form of equation (3.1), and referred to, see Ioos and Joseph (1980) as a vector fieldor ordinarydifferential equation, and to (3.1) as amap or difference equation. Both are termed dynamical systems. By a solution to Eq. (3.1) we mean amap, x, from some interval ℑ⊂ℜ1 into ℜn, which we represent as follows x : ℑ→ℜ1 t 7→ x(t) (3.3) such that x(t) satisfies (3.1), i.e. ẋ(t)= f(x(t), t;µ) (3.4) 836 K.Sibilski The dynamical systems theory (DST) provides a methodology for studying systems of ordinary differential equations. The most important ideas of DST used in the paper will be introduced in the following sections. More informa- tion onDST can be found in the book ofWiggins (1990). The first step in the DST approach is to calculate steady states of the system and their stability. The steady states can be found by setting all time derivatives equal to zero and solving the resulting set of algebraic equations. The Hartman-Grobman theorem (see p. 234 in theWiggins’s book) proves that the local stability of a steady state can be determined by linearizing the equations of motion about the steady state and calculating the eigenvalues.The implicit function theorem (Ioos and Joseph, 1980, Chap.2) proves that any steady state of a system is a continuous function of the parameters of the system if the linearized sys- tem is non-singular. A singular linearized system is characterised by a zero eigenvalue. Thus, the steady states of the equations of motion for an aircraft are continuous functions of the control surface deflections and/or vector of the thrust inclinations. Stability changes can occur as the parameters of the system are varied in such a way that the real parts of one ormore eigenvalues of the linearized system change their sign. The changes in the stability of a steady state lead to qualitatively different responses of the system and are cal- led bifurcations. Stability boundaries can be determined by searching for the steady states, which have one or more eigenvalues with zero real parts. There are many types of bifurcations, and each has different effects on the aircraft response. Qualitative changes in the response of the aircraft can be predicted by determining how many and what types of eigenvalues have zero real parts at the bifurcation point. The bifurcations for which one real eigenvalue is zero lead to the creation or destruction of two or more steady states. The bifurca- tions for which one pair of complex eigenvalues has zero real parts can lead to the creation or destruction of periodic motion. The bifurcations for which more than one real eigenvalue or more than one pair of complex eigenvalu- es has zero real parts lead to very complicated behaviour. The continuation methods are a class of a numerical algorithm used to follow a path of steady states in continuous or discrete dynamical systems as a parameter varies.They make use of the Implicit FunctionTheorem,which essentially states that if the Jacobian matrix J, Eq. (3.5), of the system linearized at a stationary point is non-singular then this solution is locally unique, i.e. it is part of a unique curve of stationary points which is a continuous function of the parameters. The parametric continuationmethods are used in numerical application of the bifurcation theory. The associated theorems involve properties of the eigenva- lues at steady state solutions points (or Floquet multipliers for periodic orbit A study of the flight dynamics helicopter... 837 solutions), and it is therefore useful in the bifurcation analysis to find all re- levant solution branches within a state-parameter space whilst evaluating the eigensystem as the algorithm proceeds. It is this characteristic of the conti- nuation methods that make them suitable for the ”global” control law design task at hand: the steady states provide a substantial amount of information about the mechanics governing the system response – including that of the closed-loop controlled system. The Jacobian matrix of an equilibrium point x0 of the vector field or the fixed point x0 is the matrix J= Df(x0)=        ∂f1 ∂x1 · ∂f1 ∂x1 ... ... ... ∂fn ∂x1 · ∂fn ∂xn        (3.5) The eigenvalues of theJacobianmatrix are important for the stability analysis. The following notations are used: Vector Fields: n0 is the number of eigenvalues of Df0 with zero real parts; n+ –numberof eigenvalues of Df0withpositive real parts; n− –number of eigenvalues of Df0 with negative real parts Maps: n0 in the number of eigenvalues of DF0 onunit circle; n+ – number of eigenvalues of DF0 outside the unit circle; n− – number of eigenvalues of DF0 inside the unit circle. An equilibrium or fixed point is called hyperbolic if n0 =0, that is, it has no overvalues on the imaginary axis. Thecontinuationmethod is used to explore thenature of the systemsteady states as a parameter varies. In this context, the steady states may refer to as standard equilibria (stationary trim points such as steady level flight, steady climbs and descents, steady turns and spins) or to periodic orbits (limit cycles such as wing rock and oscillatory spin/autorotation). The evolution of bran- ches of equilibria are computed by selecting one of the m controls/parameters as the ”free” parameter (or continuation parameter), and then solving for ẋ= f(x,µ)= 0 x∈ℜn µ ∈ℜ (3.6) where µ is one of the members of ν. In the work presented here, the continu- ation method based on that of Keller (1977) is used. The results are plotted as ”one-parameter bifurcation diagrams” – two-dimensional plots of each sta- te variable versus µ, with the line-type indicating relevant information abo- ut eigenvalue locations. Although dealing solely with the steady states, this 838 K.Sibilski information establishes a powerful insight into the mechanics governing the transient responses of the system.When designing the controllers for amano- euvrable aircraft, the graphical depiction of the underlying dynamics and how this changes in the presence of the controller, is particularly advantageous. 3.1. Poincaré map In the classical approach to differential equations, the stability of perio- dic solutions is discussed in terms of the characteristic or Floquet multipliers. The Poincaré map is a more geometrical view, which is in essence equivalent with the classical ideas. To determine the Poincaré map we need three no- tions: periodic orbit, flow, and vector field. For the purposes of this paper, it is generally sufficient to regard a differential equation as a system given by equation (3.1). We say that the vector field f generates a flow Φt : U →ℜ n, where φt(x) = φ(x,t) is a smooth function defined for all x in U and t in some interval ℑ=(a,b)⊆ℜ, and Φ satisfies Eq. (3.1) in the sense that d dt ( φ(x,t) ) t=τ = f ( φ(x,τ) ) (3.7) for all x ∈ U and τ ∈ℑ. Systems of the form (1.1), in which the vector field does not contain time explicitly, are called autonomous. Often we are given an initial condition x(0)= x0 ∈ U (3.8) in which case we seek a solution φ(x0, t) such that φ(x0, t)= x0 (3.9) In this case φ(x0, t) : ℑ → ℜ n defines a solution curve, trajectory, or orbit of differential equation (2.1) based on x0. Let γ be a periodic orbit of some flow Φt in F n arising from a nonlinear vector field f(x).We take a local cross section Ξ ⊂ℜn, of dimension n−1. Denote the point where γ intersects S by p. ThePoincaré map P : U → Σ is defined for a point q ∈ U by P(q)= φτ(q) (3.10) where τ = τ(q) is the time taken for the orbit φτ(q) based on q to the first return to S. A study of the flight dynamics helicopter... 839 3.2. Bifurcation theory For steady states of aircraft motion, very interesting phenomena appear even if onenegative real eigenvalue crosses the imaginary axiswhen the control vector varies. Two cases can be considered. • The steady state is regular, i.e.when the implicit function theoremworks and the equilibrium curve goes through a limit point. It should be no- ted that the limit point is structurally stable under uncertainties of the differential system studied. • The steady state is singular. Several equilibrium curves cross a pitchfork bifurcation point, and the bifurcation point is structurally unstable. If a pair of complex eigenvalues cross the imaginary axis, when the con- trol vector varies, Hopf’s bifurcation appears (Marsden andMcCracken, 1976; Keller, 1977; Hassard et al., 1981). Hopf’s bifurcation is another interesting bifurcation point. After crossing this point, a periodic orbit appears. Depen- ding of the nature of nonlinearities, this bifurcation may be sub-critical or supercritical. In the first case, the stable periodic orbit appears (even for large changes of the control vector). In the second case the amplitude of the orbit grows in portion to the changes of the control vector. Other domain of interest concerns the behaviour of the system when pe- riodic orbits loose their stability. Three possibilities can be concerned in this case. Firstly, a real eigenvalue crosses the point +1. Limit points are observed in this case. Secondly, a real eigenvalue crosses the point−1. Period doubling bifurcation appears in this case. In the vicinity of this point, the stable perio- dic orbit of the period T becomes unstable, and a new stable periodic orbit of the period 2T appears. This type of stability loss conducts to chaotic mo- tion. Thirdly, two conjugate eigenvalues leave the unit circle. Motion lines on a stable or unstable torus surround the unstable orbit after such a bifurcation. 3.3. Continuation technique The continuation methods are a direct result of the implicit function the- orem, which proves that steady states of a system are continuous functions of the parameters of the system except for steady states at which the linearized system is singular. The general technique is to fix all parameters except for one, and trace the steady states of the system as a function of this parameter. If one steady state of the system is known, a new steady state can be appro- ximated by linear extrapolation from the known steady state, see Doedel and Kernevez (1986). The slope of the curve at the steady state can be determined 840 K.Sibilski by taking the derivative of the equation given by setting all time derivatives equal to zero. If two steady states are known, a new steady state can be ap- proximated by linear extrapolation through the two known steady states. The stability of each steady state can be determined by calculating the eigenvalues of the linearized system. Any changes in the stability from one steady state to the next will signify a bifurcation. 3.4. Methodology scheme Taking into account experience of many researches, one can formulate the following three-stepmethodology scheme (being based on the bifurcation ana- lysis and continuation technique) for the investigation of nonlinear aircraft behaviour, see Carroll andMehra (1982), Guicheteau (1990), Jahnke andCu- lick (1994). • At the first step it is supposed that all parameters are fixed. The main aim is to search for all possible equilibria and closed orbits, and to ana- lyze their local stability. This study should be as thorough as possible. The global structure of the state space (or phase portrait) can be reve- aled after determining the asymptotic stability regions for all discovered attractors (stable equilibria and closed orbits). An approximate graphic representation plays an important role in the treating of the calculated results. • At the second step the system behaviour is predicted using the informa- tion about the evolution of the portrait with the parameter variations. The knowledge about the type of encountered bifurcation and current position with respect to the stability regions of other steadymotions are helpful for the prediction of further motion of the aircraft. The rates of parameter variations are also important for such a forecast. The faster the parameters change, the more is the difference between the steady state solution and transient motion. • Finally, the numerical simulation is used for checking the obtained pre- dictions and obtaining transient characteristics of the system dynamics for large amplitude disturbances of the state variable and parameter variations. 3.5. Steady state conditions The bifurcation theory is a set of mathematical results, which aims at the analysis and explanation of unexpectedmodifications in asymptotic behaviour of non-linear differential systems when parameters are slowly varying. For a A study of the flight dynamics helicopter... 841 fixed control vector u, two types of asymptotic states are commonly encoun- tered. The following relation gives the first one f(x,µ)= 0 (3.11) This relation is named a steady state. The second relation is given by the equation x(T)= x(0)+ T ∫ 0 f(x,µ) dt (3.12) Starting with an approximation of a steady state for given values of parame- ters, the computer code determines, by a continuation process, the solution curve x(µ) of the following set of non-linear algebraic equations, anddetermine the type of bifurcation                  Equilibriumpoints : f(x,µ)= 0 Limitpoints : f(x,µ)= 0 λ =0 Hopfpoints : f(x,µ)= 0 λ1,2 =±2iπ/T Periodicorbits : x(T)= x(0)+ T ∫ 0 f(x,µ) dt (3.13) The continuation process assumes that all functions for (3.11) are continuous and have derivatives. There are several continuation method algorithms. In the present work, the algorithms developed by Doedel and Kernevez (1986), which are based on the work of Keller (1977) were used. For the examination purposes of this paper the continuation and bifur- cation software XPPAUT (Ermentrout, 2001), aWINDOWS R© version of the well known AUTO971 software, has been used. 4. Results The helicopter configuration selected for this study is a representative of the PZL W-3 ”Sokół”. The weight without an external load is 5000kg, cor- responding to CT/σ =0.065. The altitude for all results is 0m. Most results refer to combinations of cable lengths within the range 0 < l < 55m, and 1A very useful freeware is available on the Internet: ftp://ftp.cs.concordia.ca/pub/doedel/auto. It enables one to determine all desired bifurcation points for different values of dynamical system parameters. 842 K.Sibilski load masses within 0 < m < 1500kg. The external load is assumed as a 2m×2m×5m container. Fig. 3. Steady states for the helicopter carrying the suspended load The presented results correspond to a bare airframe configuration. Figures 3 and 4 show the steady states for the helicopter with a slung load as a function of the cable length and suspended loadmass for two velocities of the system. The first one is the hovering flight (V = 0m/s), and the second one is the steady level flight (V =25m/s). The figures show thatmultiple steady states exist for most cable lengths and suspended load masses. For example, for the hovering flight case, the vertical line representing 15m of the cable length intersects four steady states. Two of them are stable the others are unstable, so the helicopter could exhibit any of these four steady states. The segment of unstable steady states (for the hovering flight case), contains the trim conditions for the cable length fromwithin 20m < l < 32mand the slung load mass 850kg < mL < 1500kg. It occurs, due to six saddle-node or Hopf’s A study of the flight dynamics helicopter... 843 Fig. 4. Steady states for the helicopter carrying the suspended load bifurcations, for the cable length of 20, 20.5, 24, 25, 28 and 38 meters. These cable lengths and slung load masses from the above mentioned ranges can be regarded as unsafe and dangerous slung load configurations. Figures 5-13 show the time simulation of motion of the helicopter- suspended load system inwhich the cable length is assumed l =20m, the su- spended loadmass mL =1050kg (for thehoveringflight case), and l =32.3m, mL = 800kg (for the level flight case). The values of parameters assumed in that simulation put the helicopterwith the slung load in the region of unstable steady states. The figures show rapidly developing aircraft oscillations (for the hovering flight case). The slung load oscillations grow slower. While conside- ring the level flight case, the observed rotorcraft and slung load motions are rather stable. The results indicate that the external load affects the helicop- ter motion. Substantial dynamic coupling occurs with the Dutch roll mode (especially in the hovering flight case). The load mode consists of the lateral and longitudinal motion of the load. The effect of the load on the phugoid is 844 K.Sibilski Fig. 5. Longitudinal deflection of the slung load Fig. 6. Lateral deflection of the slung load Fig. 7. Pitch angle of the airframe A study of the flight dynamics helicopter... 845 Fig. 8. Roll angle of the airframe Fig. 9. Longitudinal rate of the slung load Fig. 10. Lateral deflrction rate of the slung load 846 K.Sibilski Fig. 11. Pitch rate of the airframe Fig. 12. Roll rate of the airframe Fig. 13. Flapping rate of the main rotor blade A study of the flight dynamics helicopter... 847 rather small. The magnitude and frequency of those oscillations are irregular and have a chaotic character. Figures 14-18 show the Poincaré maps of selected state parameters. It can be stated that taking into consideration theunsteady rotor-blade aerodynamic model andhysteresis of aerodynamic coefficients, one encounters significant ir- regularities in the solution to equations of motion that are characteristic for chaoticmotion.When the condition for the onset of chaoticmotion is satisfied, flapping, pitching and rollingmotions appear to have chaotic oscillations. The results obtained by Tang and Dowell (1988) confirmed that unsteady aerody- namics, including deep stall phenomena together with a strongly non-linear rotorcraft model, can lead to a chaotic response of the system. Fig. 14. Longitudinal motion of the slung load (V =25m/s) – Poincarémap Fig. 15. Lateral motion of the slung load (hovering flight) – Poincaré map 848 K.Sibilski Fig. 16. Airframe phugoidmotion – Poincarémap Fig. 17. Airframe Dutch roll – Poincaré map Fig. 18. Flapping motion of the main rotor blade – Poincarémap A study of the flight dynamics helicopter... 849 5. Summary and conclusions The paper presented a study of the flight dynamics of a helicopter with an articulated rotor, carrying a suspended load. The aircraft model included the rigid body dynamics, individual flap and lag dynamics of each blade as well as the inflow dynamics. The external load was modeled as a 3-degrees of freedom pendulum suspended from a single point. The aerodynamic load was an unsteady force in the direction determined by the local airflow (defined by the angle of attack and the slip angle of the slung load). Themain aim of the study was to apply modern methods of investigation ODE for the prediction of critical configurations of the helicopter-external slung load system. Based on the carried out investigations, the following conclusions can be drawn: • The continuation and bifurcationmethods prove to be a very useful tool for analysing equations of motion of a rotorcraft carrying an external slung load. • The efficiency of the methods makes it possible to analyse complicated aerodynamic models using complete equations of motion in the entire range of the system parameters. • The knowledge of configurations of the helicopter-slung load system, which lead to bifurcation allows one to find unsafe combinations of the hanging load mass and cable length. • The need for a precise description of aerodynamic loads is the funda- mental cause of difficulties in the analysis. • Substantial dynamic coupling can occur between the Dutch roll mode and the load motion that primarily consists of lateral displacement of the load. Because of this coupling, the Dutch roll damping can decrease with a consequent deterioration of handling characteristics. References 1. Aeronautical Design StandardADS-33D, Handling Qualities Requirements for MilitaryRotorcraft, U.S. ArmyAviation andTroopCommand, St. Louis,MO, July 1994 2. ArczewskiK., PietruchaJ., 1993,MathematicalModelling of ComplexMe- chanical Systems, E. Horwood, Chichester 3. Avanzini G., DeMatteis G., 1998,Bifurcation analysis of a highly augmen- ted aircraft model, Journal of Guidance, Control & Dynamics, 20, 1 850 K.Sibilski 4. Carroll J.V.,MehraR.K., 1982,Bifurcation analysis of non-linear aircraft dynamics, Journal of Guidance Control and Dynamics, 5, 5 5. Cicolani L.S., Kanning G., Synnestvedt R., 1995, Simulation of the dy- namics of helicopter slung load systems, Journal of the American Helicopter Society, 40, 4, 44-61 6. Cicolani L.S., McCoy A.H., Tischler M.B., Tucker G.E., Gatenio P., Marmar D., 1998, Flight-time identification of a UH-60A helicopter and slung load,Proceedings of theNATORTASymposium onSystem Identification, Madrid, Spain, 10.1-10.18 7. Cliff E.M., Bailey D.B., 1975, Dynamic stability of a translating vehicle with a simple sling load, Journal of Aircraft, 12, 10, 773-777 8. Doedel E., Kernevez J.P., 1986, AUTO – Software for Continuation and Bifurcation Problems in Ordinary Differential Equations, Caltech, Pasadena 9. Dukes T.A., 1973a,Manoeuvring heavy sling loads near hover. Part I: Dam- ping the pendulousmotion, Journal of the American Helicopter Society, 18, 2, 2-11 10. Dukes T.A., 1973b,Manoeuvring heavy sling loads near hover. Part II: Some elementary manoeuvres, Journal of the American Helicopter Society, 18, 3, 17-22 11. Ermentrout B., 2001,XPPAUT5.41 – the Differential Equations Tool, Cal- tech, Pasadena 12. Etkin B., 1972,Dynamics of Atmospheric Flight, JohnWilley, N. York 13. Fusato D., Guglieri G., Celi R., 1999, Flight dynamics of an articulated rotor helicopter with an external slung load, Proceedings of the 55-th Annual Forum AHS, Montreal, Canada 14. Gabel R., Wilson G.J., 1968, Test approaches to external sling load insta- bilities, Journal of the American Helicopter Society, 13, 3, 44-55 15. Guicheteau P., 1990, Bifurcation theory in flight dynamics an application to a real combat aircraft, ICAS-90-5.10.4,Proceedings of 17th ICAS Congress, Stockholm, Sweden 16. Gutowski R., 1972,Analytical Mechanics, PWN,Warsaw 17. Hassard B.D., Kazarinoff N.D., Wan Y.H., 1981, Theory and Applica- tions of Hopf Bifurcation, Cambridge University Press 18. Ioos G., Joseph D., 1980, Elementary Stability and Bifurcation Theory, Springer-Verlag, NewYork 19. Jahnke C.C., Culick F.E.C., 1994, Application of bifurcation theory to the high-angle-of-attack dynamics of the F-14, Journal of Aircraft, 31, 1, 26-34 A study of the flight dynamics helicopter... 851 20. Keller H.B., 1977, Numerical Solution of Bifurcation and Nonlinear Eige- nvalue Problems, Application of Bifurcation Theory, Academic Press, N. York 21. Krotophalli K.R., Prasad J.V.R., Peters D.A., 1999, Helicopter rotor dynamics inflowmodelling formanoeuvring flight,Proc. 55th Annual Forum of the AHS, Montreal, Canada 22. Lucassen L.R., Sterk F.J., 1965, Dynamic stability analysis of a hovering helicopter with a sling load, Journal of the American Helicopter Society, 10, 2, 6-12 23. ŁucjanekW., SibilskiK., 1978,Longitudinaldynamic stabilityof ahovering helicopter with a hanging load, Journal of Applied and Theoretical Mechanics, 17, 2, 264-276 24. Marsden J.E., McCracken M., 1976, The Hopf bifurcation and its appli- cations,Applied Mathematical Science, 19, Springer Verlag, NewYork 25. Marusak A.J., Pietrucha J.A., Sibilski K.S., 2000, Prediction of aircraft criticalflight regimesusing continuationandbifurcationmethods,AIAAPaper, AIAA-2000-0976, 38th Aerospace Sciences Meeting and Exhibit, Reno, NV 26. Nagabhushan B.L., 1985, Low-speed stability characteristics of a helicopter with a sling load,Vertica, 9, 345-361 27. Narkiewicz J., 1994, Rotorcraft eeromechanical and aeroelastic stability, Scientific Works of Warsaw University of Technology, Issue Mechanics, 158 28. Pitt D.M., Peters D.A., 1981, Theoretical prediction of dynamic inflow derivatives,Vertica, 5, 1 29. Poli C., Cromack D., 1973, Dynamics of slung bodies using a single-point suspension system, Journal of Aircraft, 10, 2, 80-86 30. SheldonD.F., 1977,Anappreciationof thedynamicproblemsassociatedwith the external transportation of loads from a helicopter state of the art,Vertica, 1, 281-290 31. Sibilski K., 1980, Dynamic stability of a single rotor helicopter carrying a suspended payload, PhDThesis,WarsawUniversity of Technology,Warsaw 32. Sibilski K., 1998, Modelling of agile aircraft dynamics in limiting flight con- ditions, Rep. No. 2557/98,Military University of Technology,Warsaw 33. Sibilski K., 1999a, Bifurcation analysis of a helicopter non-linear dynamics, The Archive of Mechanical Engineering, 44, 2, 171-192 34. Sibilski K., 1999b, Nonlinear flight mechanics of a helicopter analysis by ap- plication of continuationmethods,Proceedings of the 25th European Rotorcraft Forum, Paper no. H-4, Roma, Italy 35. Sibilski K., 2000, An agile aircraft non-linear dynamics by continuation me- thods and bifurcation theory, ICAS-3.11.2,Proceedings of the 22nd ICAS Con- gress, Harrogate, UK 852 K.Sibilski 36. Sibilski K., Łucjanek W., 1983, Dynamic stability of a helicopter with a suspended load,The Archive of Mechanical Engineering, 30, 3-4, 249-267 37. Simpson A., Flower J.W., 1981, Lateral flutter of loads towed beneath he- licopters and its avoidance,Vertica, 5, 337-356 38. Szustak L. S., Jenney D., 1971, Control of large crane helicopters, Journal of the American Helicopter Society, 16, 3, 11-22 39. Tang D.M., Dowell E.H., 1998, Unsteady eerodynamic forces and aeroela- stic response for external store of an aircraft, J. of Aircraft, 35, 5 40. Tran C.T., Petot D., 1981, Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter rotor blade in forward flight,Vertica, 5 41. Wiggins S., 1990, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, NewYork Studium dynamiki lotu śmigłowca z podwieszonym ładunkiem z wykorzystaniem teorii bifurkacji i metod kontynuacyjnych Streszczenie W pracy przedstawiono studium dynamiki lotu śmigłowca z przegubowym wir- nikiem nośnym, przenoszącego podwieszony pod kadłubem ładunek. W zastosowa- nym modelu wiropłata uwzględniono stopnie swobody nieodkształcalnego kadłuba, dynamikę wahań i odchyleń łopat wirnika nośnego oraz dynamikę przepływu przez płaszczyznę wirnika nośnego. Założono, że podwieszony ładunek jest punktem ma- terialnym, na który działają siły aerodynamiczne, podwieszonym w jednym punkcie pod kadłubem śmigłowca. Wyniki obliczeń uzyskano dla podwieszonych ładunków o masie do 1500kg (stosunek masy ładunku do masy śmigłowca do 35%), podwie- szonych na linie o długości do 55m. Obecność podwieszonego ładunku modyfikuje charakterystyki dynamiczne i osiągi śmigłowca zewzględu na silne sprzężenia aerody- namiczne i bezwładnościowepomiędzy jego ruchema ruchami śmigłowca. Zewzględu na fakt, że układ śmigłowiec-podwieszony ładunek jest opisany za pomocą silnie nieli- niowych zwyczajnych równań różniczkowych, zastosowanie klasycznej analizymodal- nej układu nie zawsze jest możliwe. Doskonałych narzędzi do badania takich równań dostarcza teoria układów dynamicznych i będąca jej częścią teoria bifurkacji.W pra- cywykorzystanometodologię teorii układówdynamicznychdo prognozowanianatury niestabilności spowodowanej występującymi bifurkacjami. Ponadto przeprowadzono symulacje ruchu układu śmigłowiec-podwieszony ładunek po wystąpieniu bifurkacji. Manuscript received December 20, 2002; accepted for print July 14, 2003