1. Computational framework for investigation of aircraft nonlinear dynamics M.G. Goman *, A.V. Khramtsovsky De Montfort University, Leicester, England LE1 9BH, United Kingdom Received 15 October 2006; received in revised form 8 February 2007; accepted 12 February 2007 Available online 24 May 2007 Abstract A computational framework based on qualitative theory, parameter continuation and bifurcation analysis is outlined and illustrated by a number of examples for inertia coupled roll maneuvers. The focus is on the accumulation of computed results in a special database and its incorporation into the investigation process. Ways to automate the investigation of aircraft nonlinear dynamics are considered. Ó 2007 Published by Elsevier Ltd. Keywords: Nonlinear aircraft dynamics; Qualitative theory; Continuation method; Bifurcation analysis; Object oriented model; Database; Roll coupling problem 1. Introduction The increasing requirements for aircraft performance and maneuverability and the demand for more reliable pro- vision of flight safety have led to the implementation of more complex and adequate nonlinear models for modern aircraft. These models in their turn require an efficient set of computational tools for aircraft dynamics simulation and nonlinear analysis, for clearance of flight control laws with account of inherent model nonlinearities and control constraints [1]. An application of qualitative theory and bifurcation analysis methods for investigation of nonlinear aircraft dynamics is now a well established approach. Various kinds of aircraft dynamic instabilities and loss of control have been effectively investigated using continuation tech- nique and bifurcation methods. The critical flight regimes such as high incidence departures, wing rock behavior, spin and jump-like instabilities during intensive roll-coupled maneuvers have been effectively investigated within a uni- fied framework. The advanced computational methods in flight dynam- ics were first introduced in [2,3]. Several further contribu- tions presented in [4 12] proved the efficiency of the proposed methods in a number of applications for an open-loop and closed-loop aircraft dynamics. Global sta- bility and bifurcation analysis methods, continuation and eigenstructure assignment techniques have been success- fully applied to several control law design problems in [13,14]. Bifurcation diagrams for high incidence flight turned out to be a useful tool for planning of piloted sim- ulations [15,16]. A review of these research results can be found in [17,18]. Recently the bifurcation analysis method has been extended for evaluation of aircraft trim at straight-and- level flight and at level turn by adding to the problem a set of kinematical constraint equations [19,20]. This exten- sion generalizes aircraft performance evaluation using sta- bility and controllability analysis typical for one-parameter continuation techniques. A review of the latest results in the application of bifurcation and continuation methods to flight dynamics problems is given in [21]. Although considerable progress in this area has been achieved during the last two decades, there still remains a significant gap between the academic research and the acceptance of the new methods in aeronautical engineering 0965 9978/$ see front matter Ó 2007 Published by Elsevier Ltd. doi:10.1016/j.advengsoft.2007.02.004 * Corresponding author. E mail address:
[email protected] (M.G. Goman). www.elsevier.com/locate/advengsoft Available online at www.sciencedirect.com Advances in Engineering Software 39 (2008) 167 177 2. practice. The main reason for that is the lack of efficient and reliable software tools which can be adopted by users who are not experts in application of qualitative and bifur- cation analysis methods [18]. Qualitative analysis of aircraft nonlinear dynamics incorporates a number of diverse computationally-inten- sive procedures, which generate different types of data for further consideration in planning of numerical simulations, in computation of regions of attractions and in construc- tion of phase portrait representation. This paper describes the authors’ experience in developing a computational framework for the investigation of aircraft nonlinear dynamics. The framework is based on qualitative theory, a continuation method and bifurcation analysis techniques. The set of software tools outlined in this paper is an exten- sion and further development of the KRIT package [22] in a MatlabÒ environment. Some results for this computa- tional framework have been presented earlier in [23]. There was no intention in this paper to provide a user guide; it is mostly focused on some methodological aspects of the developed computational framework. A brief outline of the aircraft motion equations along with the elements of qualitative theory adopted in flight dynamics are presented in Section 2. This section also includes a description of interactive capabilities for aircraft dynamics investigation using coherent parameter continuation and phase portrait construction. The object-oriented model of the autono- mous nonlinear dynamic system used as a phase portrait database structure is outlined in Section 3. The continua- tion database features and outline of some new classes of algorithms are presented in Section 4. And finally a few computational examples for the inertia roll-coupling prob- lem are given in Section 5. 2. Goals of the qualitative investigation The motion equations suitable for bifurcation and qual- itative methods are represented in the form of an autono- mous nonlinear system of first order differential equations (the explicit form is given in the Appendix): dx dt ¼ Fðx; dÞ; x 2 X & Rn ; d 2 U & Rm ð1Þ where the state vector includes velocity vector parameters, rigid body angular rates, the pitch and bank angles x (V,a,b,p,q,r,h,/)T , the control vector includes elevator, aileron, rudder and thrust throttle deflections d (de,da,dr,dT)T . This system of equations is normally used for aircraft performance evaluation and investigation of spin dynamics. For maneuverability evaluation when aircraft perform an intensive velocity vector roll maneuver a simplified system (1) with the reduced-order state vector x (a,b,p,q,r)T and control vector d (de,da,dr)T can be considered [31,33,34]. All equilibrium states of system (1) are defined by solu- tions of the following system of nonlinear algebraic equations: Fðx; dÞ ¼ 0; x 2 X & Rn ; d 2 U & Rm ð2Þ The eigenvalues and eigenvectors of the Jacoby matrix J ¼ oF ox , calculated at equilibrium point under consideration, are used for evaluation of local stability and reconstruction of invariant manifolds of trajectories associated with equilibrium. Any equilibrium solution of system (1) with constant- control input d corresponds to steady coordinated flight along a helical trajectory with a vertical axis. Note that straight-and-level flight and the level turn belong to this class of equilibria. An aircraft in steady spin also follows a helical trajectory with small, comparable with the wing- span (R/l % 1), turn radius, with high angle of attack a % 40 80° and practically vertical flight path angle c % Àp/2. The reduced 5th-order system is used for investi- gation of the roll-coupling problem. The equilibria of this simplified system of motion equations are considered as the pseudosteady state (PSS) solutions in comparison with the full system equilibria, because the effect of the gravita- tional force is ignored and speed is assumed to be constant [31,34]. Different time-varying control inputs (1) allow one to simulate a variety of aircraft maneuvers (note that in this case system (1) should be extended by equations for posi- tion in space X, Y, Z and yaw angle w). For illustration purposes Fig. 1 shows the visualization of a simulated flight trajectory with aircraft departure at high incidence, entry into developed spin and further successful recovery from spin. Fig. 2 shows the visualization of a flight trajectory which includes intensive velocity roll. Using only direct numerical simulation the identification of maneuver limits due to onset of instability or loss of con- trol is an extremely difficult and time-expensive task. The qualitative theory, continuation and bifurcation analysis Fig. 1. High incidence departure, steady spin and spin recovery. 168 M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 3. methods facilitate this task providing estimates for critical control inputs and dangerous maneuver boundaries. Nor- mally these estimates are further verified in numerical and piloted simulation using the full mathematical model. The basic principals of the qualitative theory and bifur- cation analysis methods can be found, for example, in [24 26]. An introduction to continuation methods is given in [27]. In this paper these methods and techniques are dis- cussed in the context of their application to flight dynamics problems. At constant-control input the state space of system (1) incorporates all the trajectories representing aircraft dynamics for different initial conditions. The topology of the system state space, or the phase portrait, is defined by a set of its special trajectories, namely equilibrium points, closed orbits, toroidal manifolds and chaotic attractors. The investigation in flight dynamics is mostly focused on aircraft equilibria or trim conditions and on closed orbits representing aircraft limit cycle oscillations (LCO). A special interest is given to bifurcational parameters when the phase portrait experiences qualitative transfor- mations like changes in a number of equilibrium points, onset or disappearance of a closed orbit. Local bifurcations can be analyzed entirely via migrations of equilibrium eigenvalues and periodic orbit multipliers across the imag- inary axis or the unit circle, respectively. The most important bifurcations arising in flight dynamics applications are related to the limit point or sad- dle-node bifurcation (Fig. 3), signifying aircraft jump-like departure, and with the Hopf bifurcation or oscillatory instability, associated, for example, with the onset of wing rock at high incidence flight or flutter in aeroelastic systems (Fig. 4). Aircraft critical regimes such as spin, autorotation in roll or wing rock oscillations can coexist with the normal flight conditions at the same control inputs. Entry to one of these critical regimes depend on prehistory of control actions and also on external disturbances. This multiattrac- tor nature of aircraft dynamics requires one to perform analysis of regions of attraction (or domains of attraction) in addition to local stability for each stable regime of motion. For evaluation of regions of attraction in the pre- sented framework the method outlined in [11,13] is imple- mented. This method automates the computation and visualization of two-dimensional cross-sections of the mul- tidimensional region of attraction. Two closely interconnected tasks of qualitative analysis of aircraft nonlinear dynamics consist of: • Identification of all qualitatively different types of air- craft nonlinear behavior and specification of separating Fig. 2. Intensive velocity vector roll maneuver. Fig. 3. Equilibrium limit point or saddle node bifurcation. Fig. 4. Hopf bifurcation for equilibrium point, onset of limit cycle. M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 169 4. surfaces between corresponding regions in the parame- ter space. Analysis of their sensitivity to variation of sys- tem parameters is the important part of such an investigation. This allows the researcher to determine the safe boundaries for normal flight conditions and to isolate critical flight regimes. • Investigation of aircraft dynamics for each mode of behavior considering the phase portrait representation at different constant-control inputs and evaluation of regions of attraction for every important stable solution. The continuation procedure is the main computational tool for performing numerical bifurcation analysis [27 29,22]. It is used to study the dependence of steady-state solutions on variation of system parameters. The main fea- ture of this algorithm is the introduction of a general arc- length for the solution curve in the extended state space, which is the cartesian product of the state space and one- parameter axis (i.e. the selected bifurcation parameter). The predictor-corrector algorithm is implemented for com- putation of this solution curve. In continuation of aircraft equilibria the saddle-node bifurcation points may have a form of sharp kinks, because functions in aircraft aerodynamic models are normally made by the linear interpolation of tabulated experimental data. The continuation algorithm in this case needs special adjustment to guarantee a reliable passing through these critical points [22]. For computation of the closed orbits different versions of mapping technique are used. The presented computa- tional framework supports the Poincare´ mapping and the time-advance mappings with fixed and free periods. Vari- ous checks are performed during continuation of equilib- rium and closed orbit solutions to allow localization and classification of the encountered bifurcation points. Special procedures automatically locate, or at least try to locate, all the branches emanating from the bifurcation and branch- ing points, in order to capture the complete set of solutions. There is a number of interactive user interfaces for man- aging the investigation process using a broad spectrum of computational procedures based on qualitative methods. They include the tools for storage and retrieval of com- puted data and also special algorithms guiding the investi- gation by processing already accumulated data. Figs. 5 and 6 show the continuation and the phase por- trait GUI’s (Graphical User Interfaces), respectively. They both are used for managing the process of qualitative and bifurcation analysis. The continuation GUI includes three windows where the continuation branches in different pro- jections with local stability information can be displayed. It is possible to save the computed data in a database and to reload previously saved data. There is also a number of tools and control panels facilitating computational investi- gation, in particular, for initiation of continuation of equi- libria and closed orbits, for numerical simulation, for calculation of two-dimensional cross-sections of stability region, for opening a new or existing continuation or phase portrait database, for report generation. The phase portrait can be composed using the phase portrait GUI from the steady-state solutions available in the continuation database. This set of solutions can be sup- plemented by a number of trajectories from stable and unstable invariant manifolds associated with the equilib- rium point or the closed orbit. All this may be added by a number of regular trajectories specified for computation in the simulation panel. Fig. 5. Control panel for continuation and bifurcation analysis. 170 M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 5. 3. Object-oriented model of the autonomous nonlinear dynamic system Qualitatively different trajectories of dynamical system (1) can be considered as objects grouped in a number of classes. The object-oriented approach to system modeling provides a useful basis for development of a database struc- ture for storing computational results from nonlinear dynamic analysis. An attempt to model the autonomous nonlinear dynamical system using the class diagram adopted from the Unified Modeling Language (UML) [30] is presented in this section. The system is modelled using three classes of objects. The first class includes invariant manifolds of steady trajec- tories (IMST). It has four children classes, namely, the class of equilibrium points with dimension n 0, the class of closed orbits with dimension n 1, the class of toroidal invariant manifolds with dimension n 2,3,N À 1, where N is the order of the dynamical system, the class of strange attractors with fractional dimension p (Fig. 7). The second class includes invariant manifolds of tran- sient trajectories (IMTT) approaching or departing from the invariant manifolds of steady trajectories (IMST). For example, any invariant manifold of steady trajectories (i.e. equilibrium, closed orbit, etc.) is associated with a sta- ble invariant manifold of transient trajectories approaching it and an unstable invariant manifold of transient trajecto- ries departing from it. If IMST is stable IMTT has a full dimension of the state space and coincides with the IMST region of attraction (RA). Dimensions of a stable IMTT Fig. 6. Control panel for phase portrait investigation. Fig. 7. Object oriented model of the autonomous nonlinear dynamical system (1): the UML class diagram. M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 171 6. and an unstable IMTT coincide respectively with the num- ber of stable and unstable eigenvalues for an equilibrium and the number of stable and unstable multipliers for a closed orbit. A typical example is a separatrix surface formed by a stable manifold of trajectories approaching an aperiodically-unstable equilibrium or a closed orbit. The third class includes all regular trajectories (RT), which are different from the first and the second class of objects. They do not have association with IMST and IMTT objects, but can approach them with any required accuracy. The attributes of the objects from the first class include a state space location, eigenvalues and eigenvectors of the Jacoby matrix. The attributes of a stable object from the first class can include estimates for its region of attraction. There are two main operations with the objects. The first operation is the time advance mapping along the trajectory i.e. numerical integration of the dynamical system on some time interval, and the second operation is continuation of the object due to variation of some system parameter. 4. Database for computational results and data processing algorithms Keeping and reusing the already computed results can significantly enhance the continuation algorithm’s perfor- mance. The novelty of this approach is in the tight integra- tion between the database and the continuation algorithm. For example, at each continuation step a newly computed solution is checked against the database to determine whether such a solution has been already computed, or whether some other steady-state solution is approaching the current branch. Such a feature resulted in significant time savings, especially in the presence of complex branch- ing points or during continuation of periodic solutions. The existing databases can be used as a source of initial information for a number of algorithms. For example, a one-parameter continuation database can provide limit points for computation of two-parameter bifurcation diagram. Another example, all steady-state solutions at constant-control input can be extracted from the continua- tion database and used for phase portrait representation and if some new attractors are found using the phase por- trait investigation they can be returned back to the contin- uation database. A database structure for the one-parameter continua- tion algorithm is the most thoroughly tested one. It con- tains information about various kinds of bifurcation points and also about error points, where computation failed to converge to a steady solution, etc. Descriptive information is given for each computed and stored contin- uation branch with an outline of relationships between dif- ferent branches. The phase portrait investigation requires special compu- tational activities. Collecting data for phase portrait repre- sentation has an iterative and incremental nature and the main goal was to automate this process. A special database structure was developed for accumulating results of com- putational analysis (see Section 3) and a graphical user interface was created for managing the computational process. The phase portrait GUI, shown in Fig. 6, unifies a num- ber of different algorithms while the database is used for data exchange and storage. To create a particular phase portrait view, the following data are usually computed: • steady-state solutions and their stability (the data may be extracted from an existing continuation database or alternatively a number of various methods are provided in the phase portrait GUI for computation of steady- state solutions); • special trajectories, incoming to or outgoing from a steady-state solution along eigenvectors, corresponding to stable and unstable eigenvalues, respectively; • numerical simulation of system’s dynamics at fixed and time-varying control inputs providing time histories and visualization of aircraft spatial trajectory; • two-dimensional cross-sections of the regions of attrac- tion for all stable steady-state solutions. There is a number of auxiliary computational tools such as random search for equilibrium and periodic solutions, procedure for automatic specification of the system’s tran- sition after meeting the saddle-node bifurcation point, etc. A procedure for search of points on a surface separating two stable attractors allows finding additional unstable equilibrium and periodic solutions. A low-level compatibility between the phase portrait and the continuation databases allows a quick generation of the phase portrait representation from data already obtained by continuation. Conversely, the data collected in the phase portrait database can be exported back to the continuation database allowing continuation of new branches of solutions. 5. Computational examples A classical aircraft nonlinear problem for inertia-cou- pling roll maneuvers outlined in [31 34,21] will be consid- ered in this section to present the computational framework in action. The 5th order mathematical model for a hypothetical swept-wing fighter at high altitude super- sonic flight conditions presented and analyzed in [17] (pp. 553 559) is used for computational analysis. There is a rep- lication of some results from [17] and also some new forms such as phase portrait views, two-parameter continuation diagrams and root-loci for separate continuation branches are given and accompanied by comments on reliability and convergence of the algorithms employed, an indication of computation time and of data storage requirements. The aircraft dynamics during intensive rotation in roll is affected by inertia coupling between the longitudinal and lateral motions leading to a multitude of nonlinear phenomena such as bifurcation-induced departures, 172 M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 7. multiattractor behavior, apparent loss of controllability, hysteresis-type responses to control inputs, etc. All these phenomena can be effectively investigated using the pre- sented computational framework with a variety of tools for analysis, data storage and interactivity. The continuation of equilibrium solutions is performed as a single procedure, which automatically generates sets of initial solutions for several values of the continuation parameter and after that sequentially carries out continua- tion from these starting points. At every selected value of the continuation parameter there may be a number of start- ing points, which generate different branches of solutions. All new branches are stored in a continuation database. Special protection against recontinuation of the same branch is made by means of synchronous comparison of the current point with the data already stored in a data- base. The final set of equilibrium branches computed in an extended state space and stored in a database can be visualized in different projections. To give information about equilibrium eigenvalues the branches are drawn using different line types, for example, a solid line means that all eigenvalues are stable, a dashed line means that there is one positive real eigenvalue, a dashed-dotted line corresponds to one complex pair with positive real part, etc. All bifurcation points (i.e. the limit and the Hopf bifur- cation ones) are specially marked using information stored in a database. The obtained branches of equilibrium solutions in two projections for angle of attack vs aileron deflection and for roll rate vs aileron deflection are presented in Fig. 8. There is one branch of equilibrium solutions corresponding to normal flight conditions with small angles of attack. Two other branches correspond to aircraft autorotation regimes, which even exist at zero aileron deflection. At two limit points on these branches at da 5.5° and da À7.5°, respectively for negative and positive rotation, transitions take place to the normal branch leading to decrease in intensity of rotation (see arrows in Fig. 8). After return to the normal branch the roll rate returns to zero when aileron input is removed. At the next step the user can initiate the continuation of limit cycle oscillations (LCO) or closed orbit solutions starting from the already computed Hopf bifurcation points. The two computed families of periodical solutions are shown in Fig. 8 for positive values of aileron deflection. The amplitudes of limit cycles are superimposed on the equilibrium branches (note that similar families of period- ical solutions can be computed for negative aileron deflec- tions). The limit cycle oscillations on the main equilibrium branch become unstable after meeting the secondary Hopf bifurcation point at da 18°. Further increase of aileron input produces departure to stable limit cycle oscillations settled on the ‘‘autorotation’’ equilibrium branch with very high values of angle of attack. Based on these continuation and bifurcation analysis results the engineer can draw two main inferences. The first one is that there is an apparent loss of aileron efficiency at da > 8° (practically zero increase in angular velocity). This happens due to increasing counter-rotative aerodynamic moment produced by sideslip, which is generated by inertia coupling. The second conclusion is that at 12° < da < 18° the velocity vector roll maneuver becomes very agitated due to existence of stable limit cycles and there is also a critical aileron deflections da > 18°, when aircraft can depart to an agitated autorotation regime with very high angles of attack (see arrow in Fig. 8). Transition to high angles of attack at supersonic speeds is very dangerous as the normal and side load factors can easily exceed the structural limit. A more accurate evaluation of the aerody- namic loads during transition to high angles of attack requires inclusion into the simplified 5th-order mathemati- cal model of the nonlinear dependencies of aerodynamic characteristics and additional dynamic equation for flight velocity. For more thorough consideration of multiattractor dynamics the phase portrait views for different aileron deflections are constructed using information stored in Fig. 8. Continuation of equilibria and closed orbits: angle of attack and roll rate vs aileron deflection. M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 173 8. the continuation database. At small aileron deflections the system has five different equilibria, three stable and two aperiodically unstable. The phase portrait for da 0 is pre- sented in Fig. 9 (top plot). In addition to equilibrium points taken from the continuation database several special trajec- tories are integrated using the phase portrait GUI and added to the phase portrait. Four unstable trajectories associated with two saddle points approach three stable points. For each equilibrium point a trajectory from the associated stable invariant manifold is reconstructed by reversed time integration from its close vicinity. The stable invariant manifolds for unstable points are most important as they form the boundaries of regions of attraction for three stable equilibrium points. The phase portrait at positive aileron deflection da 14° is shown in Fig. 9 (bottom plot). There are two unstable equilibrium points and two stable solutions, a stable limit cycle on the normal equilibrium branch and a stable equi- librium regime on the autorotation branch. Unstable tra- jectories associated with saddle equilibrium point approach two stable attractors, the stable invariant mani- fold for saddle equilibrium point separates the regions of attraction for autorotation regime at high angle of attack and agitated roll maneuver at small angle of attack. Along with steady state solutions the continuation data- base stores eigenvalues and multipliers for equilibrium and closed orbits, respectively. The trajectories of eigenvalues displayed in the continuation GUI (Fig. 5) for a specified branch of solutions can be saved in a separate figure for report generation. For example, Fig. 10 shows the trajecto- ries of eigenvalues for two different branches of equilibrium solutions presented in Fig. 8. Every locally stable equilibrium and closed orbit shown in Fig. 8 has a bounded region of attraction (RA). For evaluation of the multidimensional RA a number of its two-dimensional cross-sections are generated by direct integration of dynamical system (1). Initial conditions are taken from a grid of points on the considered two-dimen- sional cross-section plane. Final destinations of integrated trajectories are identified by their entering into a close proximity of one of the available attractors stored in the Fig. 9. Phase portrait views for a, b, p projection at da = 0 (top plot) and for da = 14° (bottom plot). Fig. 10. Root loci vs aileron deflection for two branches of continuation diagram in Fig. 8. 174 M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 9. database [11,17]. For three stable equilibrium points at da 0 two cross-sections of regions of attraction in the plane (b,a) at p q r 0 and in the plane (r,p) at q a b 0 are shown in Fig. 11. Trajectories attracted by equilibrium with zero roll rate are marked by solid dia- monds. Trajectories attracted by equilibrium with positive roll rate are marked by ‘‘·’’ and trajectories attracted by equilibrium with negative roll rate are marked by ‘‘+’’. The computed cross-sections, for example, allow one to specify critical disturbances in the state space leading to departure from the normal flight conditions with zero rota- tion. A locally stable equilibrium with very small size of region of attraction should be considered as practically unstable form of motion. The saddle-node bifurcations on the continuation dia- gram (8) separate on parameter axis da segments with dif- ferent number of equilibrium points. The location of saddle-node bifurcation points in the plane of two control parameters, for example, da and dr, are defined by continu- ation curves of the extended bifurcation problem [2]: Fðx; da; dr; deÞ ¼ 0; det oFðx; da; dr; deÞ ox ¼ 0 ð3Þ The procedure for continuation of the saddle-node bifurca- tion in two-dimensional plane of (da,dr) (3) takes the initial points from the database for one-parameter continuation results. Fig. 12 shows the bifurcation diagrams in the plane of aileron and rudder deflections for two different values of elevator deflection de 0 (left plot) and de À5° (right plot). The computed results specify regions with different number of equilibria and thus allow one to assign the crit- ical boundaries in the parameter space associated with an aircraft jump-like departure. Size of the database containing the results presented in this section is about 2 Mb. The computation time for con- tinuation of three equilibrium branches is about 3 min on PC (Pentium IV, 2.4 GHz), the closed orbits require about 7 min per branch, and the computation time for two- parameter continuation diagram in Fig. 12 is about 2 min. Only the computation of two-dimensional cross-sec- tions of region of attraction with a fine grid of points remains rather time-consuming it takes about 20 min to Fig. 11. Cross sections of regions of attraction for three equilibrium points in (b,a) and (r,p) planes. Fig. 12. Two parameter bifurcation diagrams in the plane of aileron and rudder deflection for a saddle node bifurcation point: regions with different number of equilibria (left plot for de = 0, right plot for de = 5°). M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 175 10. compute one cross-section. The developed computational tools were tested during investigations of nonlinear dynam- ics for a number of realistic aircraft models covering the full flight envelope. Acceptable convergence and robustness properties of the algorithms were demonstrated, as well as the possibility to restart computations from a checkpoint in case of hardware or software failure. 6. Conclusions The presented computational framework provides a full suite of algorithms and methods required for investigation of aircraft nonlinear dynamics by a systematic analysis of aircraft equilibrium and periodical forms of motion. The computational tools for one- and two-parameter continua- tion, bifurcation analysis, phase portrait/time histories, evaluation of regions of attraction, root loci, systematic starting solution solver, etc. are linked together in a logical manner that promotes a structured approach to problems by engineers who are not necessarily highly specialized in the field of qualitative methods and continuation software. The developed databases for accumulation of computed results, their integration with computational algorithms and an interactive way of managing the investigation pro- cess via a number of Graphical User Interfaces provide a productive and user-friendly environment. The way in which computed information is accessed and processed is a significant step forward in making such software user- friendly and suited to an engineering practice. Acknowledgements The presented computational framework for aircraft nonlinear dynamics investigation was developed during the research project funded by DERA/QinetiQ Ltd., Bed- ford, UK. During various years it was monitored by Phill Smith, Darren Littleboy and Yoge Patel. Special thanks go to the reviewers and also to Daniel Walker for their valuable and constructive comments and suggestions, which helped to improve the presentation. The authors acknowledge the contribution to the computational system of Evgeny Kolesnikov. Appendix. Aircraft equations of motion The autonomous set of equations for rigid aircraft flight dynamics [2]: _V ¼ ð1=mÞ½ððX þ TÞ cos a þ Z sin aÞ þ Y sin bŠ Á Á Á þ gðsin b cos h sin h À cos a cos b sin h þ sin a cos b cos h cos /Þ _a ¼ q À ðp cos a þ r sin aÞ tan b þ ð1=mV Þ Â ½Z cos a À ðX þ TÞ sin aŠ Á Á Á þ ðg=V cos bÞðsin a sin h þ cos a cos h cos /Þ _b ¼ p sin a À r cos a þ ð1=mV Þ Â fY cos b À ½ðX þ TÞ cos a þ Z sin aŠ sin bg Á Á Á þ ðg=V Þðcos a sin b sin h þ cos b cos h sin / À sin a sin b cos h cos /Þ _p ¼ ððIYY À IZZÞ=IXX Þqr þ L=IXX _q ¼ ððIZZ À IXX Þ=IYY Þpr þ M=IYY _r ¼ ððIXX À IYY Þ=IZZÞpq þ N=IZZ _h ¼ q cos / À r sin / _/ ¼ p þ ðq sin / þ r cos /Þ tan h The aerodynamic coefficients CX,Y,Z 1 and Cl,m,n 2 are repre- sented as functions of motion parameters a, b, p, q, r and control deflections de, da, dr. Note that at high angles of at- tack the aerodynamic dependencies are nonlinear with a strong effect of prehistory of motion [35]. The reduced set of equations for analysis of aircraft intensive rolling maneuvers [34]: _a ¼ q À pb þ zaa þ zde de _b ¼ pa À r þ ybb þ ydr dr _p ¼ ððIYY À IZZÞ=IXX Þqr þ lbb þ lpp þ lrr þ lda da þ ldr dr _q ¼ ððIZZ À IXX Þ=IYY Þpr þ maa þ mqq þ ma _a þ mde de _r ¼ ððIXX À IYY Þ=IZZÞpq þ nbb þ npp þ nrr þ nda da þ ndr dr References [1] Fielding C, Varga A, Bennani S, Selier M, editors. Advanced techniques for clearance of flight control laws. Lecturer notes on control and information sciences, vol. 283. Springer Verlag; 2002. [2] Mehra RK, Kessel WC, Carroll JV. Global stability and control analysis of aircraft at high angles of attack, Annual Technical Reports 1/2/3, ONR CR215 (1/2/3), 1977/1978/1979, Scientific Sys tems Inc., USA. [3] Carroll JV, Mehra RK. Bifurcation analysis of nonlinear aircraft dynamics. J Guid Control Dynam 1982;5(5):529 36. [4] Guicheteau P. Application de la the´orie des bifurcations a l’e´tude del pertes de controle sur avion de combat. La Recherche Aerospatiale 1982;2:61 73. [5] Zagaynov GI, Goman MG. Bifurcation analysis of critical aircraft flight regimes. In: Proceedings of the 14th congress of ICAS 84 4.2.1, Toulouse, France, 1984. [6] Jahnke CC, Culick FEC. Application of bifurcation theory to the high angle of attack dynamics of the F 14. J Aircraft 1994;31(1): 26 34. [7] Lowenberg MH. Bifurcation analysis as a tool for post departure stability enhancement. AIAA 1997 3716, AIAA atmospheric flight mechanics conference, New Orlean, LA, 1997. [8] Guicheteau P. Bifurcation theory: a tool for nonlinear flight dynamics. Philos Trans Roy Soc London A 1998;356(1745):2181 201. [9] Planeaux JB, Beck JA, Baumann DD. Bifurcation analysis of a model fighter aircraft with control augmentation, AIAA Paper, 1990; 1990 2836. [10] Avanzini G, de Matteis G. Bifurcation analysis of a highly augmented aircraft model. J Guid Control Dynam 1997;20(4):754 9. 1 X qSCX ; Y qSCY ; Z qSCZ. 2 L qSbCl; M qScCm; N qSbCn. 176 M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 11. [11] Goman MG, Khramtsovsky AV. Global stability analysis of nonlin ear aircraft dynamics. AIAA Paper 97 3721, AIAA guidance, navigation, and control conference and exhibit, New Orleans, LA, 1997. [12] Littleboy DM, Smith PR. Bifurcation analysis of a high incidence aircraft with nonlinear dynamic inversion control. AIAA Paper 97 3717, AIAA atmospheric flight mechanics conference and exhibit, New Orleans, LA, 1997. [13] Goman MG, Khramtsovsky AV. Application of bifurcation and continuation methods for an aircraft control law design. the theme issue flight dynamics of high performance manoeuvrable aircraft. Philos Trans Roy Soc London A 1998;356(1745):2277 95. [14] Charles G, Lowenberg M, Stoten D, Wang X, di Bernardo M. Aircraft flight dynamics analysis and controller design using bifur cation tailoring. AIAA 2002 4751, AIAA guidance, navigation, and control conference and exhibit, Monterey, CA, 2002. [15] Patel Y, Littleboy D. Piloted simulation tools for aircraft departure analysis. Philos Trans Roy Soc London A 1998;356(1745):2203 21. [16] Lowenberg MH, Patel Y. Use of bifurcation diagrams in piloted test procedures. Aeronaut J 2000;104:225 35. [17] Goman MG, Zagaynov GI, Khramtsovsky AV. Application of bifurcation theory to nonlinear flight dynamics problems. Prog Aerosp Sci 1997;33(9):539 86. [18] Macmillen FBJ, Thompson JMT. Bifurcation analysis in the flight dynamics design process? A view from the aircraft industry. Philos Trans Roy Soc London A 1998;356(1745):2321 33. [19] Ananthkrishnan N, Sinha NK. Level flight trim and stability analysis using an extended bifurcation and continuation procedure. J Guid Control Dynam 2001;24(6):1225 8. [20] Paranjape A, Ananthkrishnan N. Airplane level turn performance, including stability constraints, using EBAC method. AIAA 2005 5898, AIAA atmospheric flight mechanics conference and exhibit, San Francisco, CA, 2005. [21] Paranjape A, Sinha NK, Ananthkrishnan N. Use of bifurcation and continuation methods for aircraft trim and stability analysis a state of the art. AIAA 2007 1051, 45th AIAA aerospace sciences meeting and exhibit, 8 11 January 2007, Reno, Nevada. [22] Goman MG, Khramtsovsky AV. KRIT: scientific package for continuation and bifurcation analysis for flight dynamics applica tions, TsAGI Contract Report for DRA, Bedford, 1993. [23] Goman MG, Patel Y, Khramtsovsky AV. Flight clearance tools using a non linear bifurcation analysis framework. AIAA Paper 2003 5557, AIAA guidance, navigation, and control conference and exhibit, Austin, TX, 2003. [24] Arnold VI. Geometrical methods in the theory of ordinary differential equations. Springer Verlag; 1998. [25] Kuznetsov YuA. Elements of applied bifurcation theory. Springer Verlag; 1998. [26] Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems, and bifurcations of vector fields. Applied Mathematical Sciences, vol. 42. Westview Press; 2002. [27] Allgower EL, Georg K. Numerical continuation methods: an introduction. Berlin Heidelberg: Springer Verlag; 1990. [28] Doedel EJ, Paffenroth RC, Champneys AR, Fairgrieve TF, Kuznet sov YA, Sandstede B, et al. AUTO2000: Continuation and bifurca tion software for ordinary differential equations (with HomCont), Technical Report, California Inst. of Technology, Pasadena, CA, USA, 2001. [29] Dhooge A, Govaerts W, Kuznetsov YuA. MATCONT: a MATLAB package for numerical bifurcation analysis of ODEs. ACM Trans Math Software 2003;29:141 64. [30] Fowler Martin, Kendall Scott. UML Distilled. A Breif Guide to the Standard Object Modeling Language, 2nd edition, Foreworded by Grady Booch, Ivar Jacobson, and James Rumkbaugh, Addison Wesley; 2000. [31] Hacker T, Oprisiu C. A discussion of the roll coupling problem. Prog Aerosp Sci 1974;15:151 81. [32] Schy AA, Hannah ME. Prediction of jump phenomena in roll coupled maneuvers of airplanes. J Aircraft 1977;14(4):375 82. [33] Hacker T. Constant control rolling maneuver. J Guid Control 1978;1(5):313 8. [34] Young JW, Schy AA, Jonson KG. Pseudosteady state analysis of nonlinear aircraft maneuvers, NASA TP 1758, 1980. [35] Abramov NB, Goman MG, Khrabrov AN. Aircraft dynamics at high incidence flight with account of unsteady aerodynamic effects. AIAA Paper 2004 5274, AIAA atmospheric flight mechanics conference and exhibit, August 2004, Providence, Rhode Island, USA. M.G. Goman, A.V. Khramtsovsky / Advances in Engineering Software 39 (2008) 167 177 177