## Abstract

Bifurcation analysis is a powerful method for studying the steady-state nonlinear dynamics of systems. Software tools exist for the numerical continuation of steady-state solutions as parameters of the system are varied. These tools make it possible to generate ‘maps of solutions’ in an efficient way that provide valuable insight into the overall dynamic behaviour of a system and potentially to influence the design process. While this approach has been employed in the military aircraft control community to understand the effectiveness of controllers, the use of bifurcation analysis in the wider aircraft industry is yet limited. This paper reports progress on how bifurcation analysis can play a role as part of the design process for passenger aircraft.

## 1. Introduction

Aeronautical vehicles are continually driven by market needs to provide highly reliable aircraft with better performance and functionality at a reduced cost of ownership. These characteristics form part of the emergent properties, capabilities and behaviours of the vehicle. They arise through the interactions among multiple components (internal influences) that constitute the vehicle, and the vehicle’s interactions with its operating environment (external influences). Engineers strive to design aircraft that, when built, will reliably exhibit these desirable behaviours. The use of optimization techniques to reduce weight, loads, cost of manufacture and cost of ownership, while increasing the availability and maintainability of the product, have resulted in an unprecedented level of design integration. The resulting solutions inevitably give rise to new interactions that are nonlinear.

Architects of aeronautical vehicles tend to use simplified models in the conceptual design stage, which comprise either surrogate models or linearized differential equation models, as depicted in figure 1. Surrogate models are used for determining the vehicle’s dimensions and do not contain any dynamics. Linearized models of the vehicle are used to perform a large number of simulation studies with different model parameters and operating scenarios. Although the mathematical models could be analytical, there is a tendency for the engineering analysts to resort to numerical methods by using modelling environments such as Matlab Simulink. The resulting time series are analysed to characterize the behaviours, with no guarantees that all possible behaviours will be discovered, due to the linear model assumption.

The introduction of nonlinear analyses, and more specifically bifurcation analysis, can provide the detailed dynamics information an architect may need during the concept phase. A timely understanding of nonlinear behaviours allows vehicle architects to specify architectures that lead to high-quality products, by either diminishing the likelihood of undesired nonlinear behaviours occurring, or by using nonlinear behaviours to their advantage. This is important for the vehicle designers and manufacturers—while only 8% of the total product budget is spent by the end of the concept stage, 80% of the cost of the product will have been committed [1].

This paper considers the use of bifurcation analysis and numerical continuation in this crucial phase of the design process. Bifurcation analysis has been used widely within the military aircraft control community for many years; see, for example, the reviews by Thompson & Macmillen [2] and Paranjape *et al.* [3]. Here, we concentrate on work conducted on passenger aircraft, both in flight and during ground manoeuvres. Section 2 gives a brief description of bifurcation analysis along with its software implementation. Section 3 discusses with a number of examples the insights that bifurcation analysis can provide. A more detailed industrial case study is then considered in §4, and an outlook and conclusion are given is §5.

## 2. Bifurcation analysis

Mathematical models that one encounters in the context of aerospace applications can often be written as a system of first-order differential equations of the form
2.1Here, *x*∈*R*^{n} is the vector of *n* states that describe the system, λ∈*R*^{m} is a vector of all system parameters, and the function *f* is assumed to be sufficiently smooth. System (2.1) is also referred to as a vector field, and its dynamics are given by the corresponding flow *ϕ*^{t} that maps an initial condition *x*_{0} to *ϕ*^{t}(*x*_{0},λ) after time *t*. The basic issue is to find out how the observed dynamics depend on the system parameters. As the function *f* is nonlinear in general, system (2.1) cannot be solved analytically. Hence, numerical methods need to be employed. One option is to use simulation by numerical integration from some chosen initial condition. However, it is often computationally expensive to find all behaviour of interest in this way.

An alternative is offered by the technique of bifurcation analysis. Its underlying idea is to follow or continue steady states and periodic solutions of (2.1) in a chosen continuation parameter. Recall that a steady state or equilibrium is a point *x*_{0}, where *f*(*x*_{0},λ)=0, while a periodic orbit is defined by an initial condition *x*_{0} that returns back to itself after the period *T*, that is, *ϕ*^{T}(*x*_{0},λ)=*x*_{0}. The stability of an equilibrium is determined by the eigenvalues of its Jacobian matrix; similarly, the stability of a periodic orbit is given by its Floquet multipliers (which are the eigenvalues of the return map along the periodic orbit). Hence, stability properties of equilibria and periodic orbits can be monitored during a continuation. A change of stability leads to a qualitative change of the dynamics—a bifurcation. Bifurcation points can be detected during the continuation of an equilibrium or periodic solution. From such points of bifurcations, other solutions may emerge or bifurcate, and they can then also be followed (irrespective of whether they are stable or not) to obtain a one-parameter bifurcation diagram that represents all solutions and their stability as a function of the chosen continuation parameter. Such a bifurcation diagram can be represented by projection on any of the *n* states in the vector *x*. Detected bifurcation points can be continued in an additional parameter to yield bifurcation curves in a parameter plane.

Bifurcation analysis is based on two important ingredients:

—

*Bifurcation theory*. Bifurcations can be identified and listed in a systematic way, and they are represented by normal forms, meaning that the nearby dynamics are known*a priori*.—

*Numerical continuation*. Continuation of equilibria and periodic orbits, as well as the detection and continuation of their stability changes, can be performed with available software packages; for example, with AUTO [4,5], MatCont [6] or CoCo [7].

Bifurcation analysis is a powerful tool as it allows one to identify in a systematic way where dynamics of interest exist in parameter space. Regions of stable solutions can be determined by computing stability boundaries directly as functions of the relevant system parameters. In particular, regions of multi-stability can be identified in this way. In this paper, we will encounter the following bifurcations of equilibria and periodic orbits, which are typical in the sense that they must be expected to occur when a single parameter is changed; see, for example, [8–10] for more details.

— At a

*saddle node bifurcation*, two equilibria come together and disappear; this is characterized by a single eigenvalue 0.— At a

*Hopf bifurcation*, an equilibrium loses stability and a periodic orbit emerges; this is characterized by a pair of purely imaginary eigenvalues ±i*ω*.— At a

*saddle–node of limit cycle bifurcation*, two periodic orbits come together and disappear; this is characterized by a single Floquet multiplier 1.— At a

*period-doubling bifurcation*, a periodic orbit loses its stability and a periodic orbit with twice the period emerges; this is characterized by a single Floquet multiplier −1.— At a

*torus bifurcation*, a periodic orbit loses stability and an invariant torus emerges; this is characterized by a pair of imaginary Floquet multipliers e^{±i2πω}.

The practical use of bifurcation theory requires the knowledge of the basic bifurcations that must be expected to occur, as well as familiarity with the governing physical phenomena and their modelling in software. Moreover, the model under consideration must satisfy certain minimal conditions for this approach to work. In practical terms, the right-hand side (the function *f*) must depend smoothly on the state *x* and the parameter λ; for example, data obtained from look-up tables to define *f* must be interpolated smoothly.

### (a) Software packages and the Dynamical Systems Toolbox

The continuation of equilibria, periodic orbits and their bifurcations is a well-established technique that has been implemented in a number of software packages; see, for example, [11] as an entry point to the extensive literature on the subject. Arguably, the most extensively used package is AUTO [4,5]; it was developed by Doedel in the early 1980s and has been extended over the years with additional capabilities. More recent developments are MatCont [6], which is written in Matlab and also implements normal form calculations, and CoCo [7], which facilitates the formulation of multi-segment boundary value problems.

The efficiency of the implemented continuation routines depends on the size of the system, that is, on the number *n* of states represented in the state vector *x*, as well as on how costly it is to evaluate the right-hand side *f*. Bifurcation analysis is ideally suited to reasonably small systems, say, with up to 20 states, where the function *f* is given in explicit functional form. It is more challenging but also possible to employ bifurcation analysis to larger or industry developed models.

To facilitate the use of continuation methods in an industrial context, the Development of the Dynamical Systems Toolbox (DST) was begun in 2008 at the University of Bristol, UK. Building on previous work by Ryan Bedford at the university, the goal was to incorporate the continuation software AUTO into the user-friendly environment of Matlab [12] to allow for the convenient bifurcation analysis of industrial Simulink and SimMechanics models. In this way, industrially validated models can be investigated directly without the need to convert (and often simplify) models to a format that can be used by stand-alone AUTO. In addition, using such validated models provides confidence in the results—an important part in any certification process.

The DST integrates the Fortran AUTO code into Matlab via mex-functions, ensuring usability, computational speed, toolbox support, etc. [12]. In effect, the speed of Fortran is combined with the user-friendly interface of Matlab in such a way that AUTO has direct access to the states of Simulink or SimMechanics models. Similar output objects to that of AUTO output files are generated, and any number of additional outputs can be obtained from Simulink output ports, which are written alongside the continuation parameters and states. Figure 2 depicts the architecture of the DST under Matlab.

During a recent project, LMS Virtual.Lab Motion was linked to the DST for the analysis of nose landing gear, while current work is focusing on the generation of C-code from a Simulink model in a format that can be linked in to the DST. This will avoid the use of the relatively slow `mexCallMatlab` command, when compared with pure Fortran or C-code. A C-code function of the model also allows one to split the calculations across different processors, supporting a current feature of AUTO, and would increase the calculation speed of periodic orbit continuation considerably.

More widespread use of the DST is promoted by providing documentation and reference material that is easy to use, while concrete examples act as additional training material for the user. We have combined most of the user manual of AUTO into the toolbox, which is integrated into the Matlab help environment. The DST therefore feels like any other toolbox that has been developed for Matlab; it can be downloaded from Coetzee *et al.* [13].

## 3. Applications of bifurcation analysis

We now discuss the usefulness of bifurcation analysis in terms of behaviour characterization, determining sudden changes in behaviour and supporting conceptual design. To this end, we consider some specific examples of models arising from industrial practice.

### (a) Behaviour characterization

Nonlinear dynamics result in various types of steady-state responses. It has been shown by considering nonlinear normal modes that—in stark contrast to a linear system—a two-degrees-of-freedom nonlinear oscillator can exhibit more than two modes of vibration (e.g. Kerschen *et al.* [14] or Neild *et al.* [15]).

Bifurcation analysis is the tool of choice for identifying in a convenient way the types of behaviour exhibited by a system as a parameter varies. The ability of numerical continuation to follow unstable branches, as well as stable ones, is particularly useful for identifying regions of stable solution and their boundaries. This also allows one to find regions of multiple stable solutions, without the need to run simulations for many initial conditions or for different applied perturbations to trigger a jump to another solution branch. We now demonstrate this capability with the example of *landing gear shimmy oscillations*.

Its known that pilots report cockpit vibrations under certain conditions, which could be attributed to several factors, one of which is nose-gear vibration, better known as shimmy. These vibrations are often caused by a combination of factors, which include tyre flat-spots, pressure differences between tyres, faulty shimmy dampers and excessive backlash in the gear assembly. The root cause is not always easy to determine and, therefore, a list of maintenance actions are applied by the aircraft operator. Such events have been linked to gear failure; see, for example, the accident report by CIAIAC [16].

We consider here shimmy dynamics of landing gear during forward motion. Shimmy has been studied since the 1920s with the key role of flexibility, from either the pneumatic tyres [17] or structural flexibility [18], in the onset of automobile shimmy identified. In aircraft, initially linear techniques were used to identify the onset of shimmy in nose landing gear [19]. Somieski [20] introduced nonlinearity and identified a Hopf bifurcation to stable limit cycle shimmy oscillations. The sources of nonlinearities are primarily the tyre dynamics, which are challenging to model and often add additional states and even delays to the system [21], and geometric coupling. Main landing gear, which is less prone to shimmy, has also been studied and Van der Valk & Pacejka [22] presented a case study into a main gear failure.

The typical low-order lumped parameter model of landing gear includes two degrees of freedom: the torsional rotation *ψ* and the lateral bending angle *δ*. When the fuselage is modelled as a vertical load and an additional lateral degree-of-freedom, *y*, and the tyre dynamics are described by the stretched string approach [23], the model becomes a first-order differential equation of the form (2.1) with *n*=7 states. Figure 3 shows a one-parameter bifurcation diagram of this model in forward velocity *V* , displayed in terms of torsional rotation *ψ* and lateral bending *δ*; see Terkovics *et al.* [24] for full details. In this and further bifurcation diagrams in this paper, the solution types and bifurcations are represented as shown in table 1.

Figure 3 shows that, at low and at high speed, the straightforward rolling equilibrium solution, for which there are no torsional or lateral deflections, is stable. As the velocity *V* increases from a low level or reduces from a high level, this desirable solution becomes unstable at Hopf bifurcations (at *V* ≈4.5 m s^{−1} and *V* ≈180.0 m s^{−1}, respectively) that lead to stable periodic oscillations of the landing gear. From the Hopf bifurcation point at *V* ≈4.5 m s^{−1}, a branch of single-frequency periodic orbits emerges that features a significant torsional component with little lateral motion—referred to as *torsional shimmy* [25]. This branch ends on the unstable no-shimmy solution at *V* ≈75.6 m s^{−1}. The Hopf bifurcation point at *V* ≈180.0 m s^{−1} that is encountered as *V* is reduced from a high level also leads to a branch of periodic orbits that is now dominated by lateral motion—we speak of *lateral shimmy*. This branch ends on the unstable no-shimmy solution at *V* ≈6.5 m s^{−1}.

On both of these branches, we find points of torus bifurcations where the respective single-frequency periodic orbit loses its stability. From these torus bifurcation points two branches of stable multiple-frequency periodic solutions have been identified in which both lateral and torsional motion occurs; these were calculated via successive time integrations, as tori cannot readily be computed by continuation. In addition, there exists an isola, part of which is visible in the panels of figure 3, which has a stable region spanning approximately 31.2 m s^{−1} <*V* <107.5 m s^{−1} between a saddle–node and a torus bifurcation. This stable periodic solution is characterized by significant fuselage lateral motion in addition to landing gear deflections.

Overall, a one-parameter bifurcation diagram such as figure 3 is very useful for identifying both the possible types of observed behaviour as well as the transitions between them. In this study, the parameters were selected for illustrative purposes and do not reflect real behaviour. Rather, the goal was to investigate qualitative features of nose landing gear shimmy to aid the design process, which ultimately must ensure that the onset of shimmy is well outside the operational envelope.

### (b) Sudden changes in behaviour

There are three main ways in which a steady-state response of a system can change. One is a perturbation in the external loading, for example a gust loading during steady flight. A second is a change in one of the operating parameters such as the forward aircraft velocity during a ground manoeuvre. The third is a change in a parameter of the structure; this may be slow, such as increasing freeplay of a control surface due to wear, or fast, such as stiffness changes in a tyre due to temperature increase on spin-up during landing.

Understanding how a system will respond to loading or parameter changes is often an important part of the design process. It requires knowledge of the steady-state solutions, of the likelihood of being perturbed onto a different solution branch, and of the transient response. One-parameter bifurcation diagrams provide some basic important information in this regard, and they can be used to identify parameter regions of interest that can then be studied by means of time-integration techniques to find the transient response.

To illustrate these ideas, we consider the *flight dynamics of a civil aircraft* that encounters ‘upset conditions’. Recent experience has shown that a systematic evaluation of the nonlinear dynamics can be crucial to appreciating the causes and consequences of upset. The phenomenon of airliner upset and loss-of-control (LOC) is the major cause of fatalities in commercial aviation [26]. LOC may be defined as a state outside the normal operating flight envelope that cannot be changed in a predictable manner by pilot inputs, is strongly nonlinear (e.g. due to aerodynamic load characteristics at high angles of incidence, or inertial coupling), and is likely to result in high angular displacements and rates [27]. LOC is typically triggered by faults, external events and/or inappropriate pilot inputs; it has occurred both on older aircraft without envelope-protecting control laws as well as on newer fly-by-wire types.

Mathematical models developed by airliner manufacturers normally only cover a limited angle-of-incidence range, rendering them incapable of representing behaviour in upset/LOC. However, as part of NASA’s Aviation Safety Program, a wide envelope mathematical model of a sub-scale representative twin-underwing engine airliner was created by using both wind tunnel models and the remotely piloted ‘AirSTAR’ vehicle [28,29]. This simulation tool of the sub-scale aircraft is known the Generic Transport Model (GTM), and it has been used extensively in upset/LOC studies. The recent application of bifurcation analysis to the open-loop GTM [30,31] has provided an enhanced understanding of its upset/LOC behaviour: the resulting bifurcation diagrams are essentially a ‘map’ of the vehicle’s various nonlinear responses.

Figure 4 shows an example of a one-parameter bifurcation diagram of the GTM, where the continuation parameter is elevator angle *δ*_{e}, with all other parameters fixed, and the solution is represented by the bank angle state *ϕ*. Note that decreasing elevator corresponds to the pilot pulling back on the stick, thus increasing the angle of attack, *α*. The bifurcation diagram shows steady-state and periodic solutions and their stability, allowing one to distinguish the ranges of *δ*_{e} that correspond to different observed behaviour. In region A for *δ*_{e}>−1°, one finds stable trimmed flight, where the bank angle remains zero throughout the elevator range in the absence of lateral–directional control surface inputs. Below this value of elevator *δ*_{e}, the solutions at *ϕ*≈0 are unstable (in the ‘slow spiral’ mode), all the way to *δ*_{e}>−30°, and one finds additional solution branches that correspond to different dynamics as listed in the caption of figure 4. Of those, most prominent are the two steady-state branches marked C^{+} and C^{−}, designating positive and negative bank angle solutions, respectively.

A pilot whose aircraft has encountered a spin is likely to command a wings-level condition in order to recover normal flight. Figure 4 suggests that it is possible to recover the aircraft from any of these undesirable situations by pushing the stick forward so that the elevator is in the normal flight operating region A. Figure 5 illustrates this recovery action, with a time-stepping simulation, in the plane of sideslip angle *β* (which is zero or small in conventional flight) and *α*, and the resulting aircraft trajectory. The model is initiated in the period-one orbit that constitutes oscillatory spin F, with *δ*_{e}=−24°; here, *α* oscillates between 33° and 38° and *β* between 2° and 10°. After 10 s, the elevator is ramped quickly to *δ*_{e}=2.58° in the normal trimmed flight regime A. The aircraft transitions very rapidly to this regime (*α*≈3°, *β*≈0°), although it then undergoes a low-frequency oscillation as it settles on the desired solution. The loss of altitude arising from the steep oscillatory spin is evident in the trajectory plot; this descent is then arrested by the recovery action. Hence, the recovery action suggested by the one-parameter bifurcation diagram has been confirmed in this case: the transition to recovery is indeed possible and does not endanger the aircraft further.

### (c) Supporting conceptual design

Bifurcation analysis can be used to give insight into the dominant parameters that govern the response. In this way, it can give weight to both the justification of modelling assumptions (such as the conditions under which fuselage motion must be included in a landing gear shimmy analysis [24]) and lead to the identification of operating regions or parameters that give a desirable response. This process often involves continuing bifurcation points in two parameters to obtain two-parameter bifurcation diagrams.

As an illustration, we consider the *locking mechanism of a main landing gear*. The gear consists of a main shock strut which is supported by one (or two) folding sidestays (figure 6). The sidestay is locked in the deployed state by lock-links which hit stops when the over-centre angle *θ* exceeds a certain positive value. While there is much literature on the dynamics of landing gear as a structure during landing [32], there is little research on the gear as a mechanism during deployment and this is limited to the gear kinematics [33,34]. However, the quasi-static motion of the mechanism can be formulated as a differential algebraic equation and then transformed into ordinary differential equations though differentiation of the constraints [35]. Bifurcation analysis of these equations allows the investigation of how the lock-link force varies with over-centre angle *θ*. A fold in this curve leads to a jump in over-centre angle as it ‘snaps’ through 0° and hits the stops, hence locking the mechanism [36]—providing an example of where a sudden jump is advantageous. To understand the robustness of the mechanism, this fold point can then be continued in other parameters, such as attachment point deflections due to wing flexure and compliance in the side stays [37]. Such an investigation is able to provide valuable design requirements during the conceptual stage when the design has not been frozen.

## 4. Industrial case study: aircraft ground dynamics

Nonlinear behaviour associated with landing gears can be observed in wheels during braking, as stiction in the shock absorbers and, most importantly, in the tyres during ground manoeuvres. Experience has shown that the use of different tyres can significantly change the handling qualities of an aircraft on the ground, which is very similar to what is seen on Formula One cars. Comprehensive nonlinear simulations are therefore used to ensure that an aircraft with good handling qualities is delivered to the final customer. This does however come at a high computational cost. Complementary methods from bifurcation theory overcome this problem due to the efficiency of the method, which provides an overall map of the possible dynamics. Points of specific interest can be identified for further investigation by simulations.

The starting point of our investigation of ground handling was an MSC.ADAMS model that was developed by the Landing Gear Group of Airbus for ground-handling studies [38]. MSC.ADAMS is a software program that can analyse kinematic, quasi-static and dynamic mechanical systems. The first modelling step is to describe the rigid parts and the joints connecting the parts [39], where a part is described by its mass, inertia and orientation. The nose gear is constrained by a cylindrical joint, driven by an angular motion, and the main gears are constrained by translational joints. The next step is the addition of internal force elements, known as line of sight forces, to represent the shock absorbers and tyre forces. The tyres are modelled with impact functions that switch on as soon as the distance between the wheel centre and the tyre becomes less than the wheel radius, where the lateral force on a tyre is dependent on the slip angle and the vertical force on the tyre. External forces such as thrust and aerodynamic forces, known as action-only forces, are then added. The value of such a software tool for a company lies in the fact that the equations of motion are constructed automatically, allowing the engineer to focus on the engineering aspects of the system, and not the derivation of complex equations. As MSC.ADAMS cannot be coupled to the DST, a process was developed to translate the MSC.ADAMS model to SimMechanics, which can be coupled to the DST. While there are important differences between the two modelling packages, the development of the SimMechanics model followed that of the ADAMS model.

After the successful coupling of the industrial ground-handling model to the DST, the first step in the analysis is to find an equilibrium condition for the aircraft. In our case, a velocity controller is used to accelerate the aircraft to a desired steady-state velocity over a few seconds. The thrust is then held constant until the steady-state response is reached. The corresponding values are then used as the starting conditions for the continuation analysis, which follows the equilibria as the steering angle *δ* is varied while detecting bifurcations. In this way, the dynamics of the aircraft can be investigated when it attempts a turn at a given steering angle.

Figure 7 shows a one-parameter bifurcation diagram in *δ* and corresponding aircraft ground tracks that illustrate the overall complexity of possible ground dynamics of a single-aisle aircraft; this diagram is far from the normal operating envelope of the aircraft, but does provide interesting insights into the possible observed behaviours and transitions between them due to a high thrust setting. In the one-parameter bifurcation diagram, there are four separate stable parts of equilibrium branches, as well as a branch of oscillations that is bounded by two points of Hopf bifurcation (indicated by circles in the figure).

The stable section (*a*) of the equilibrium branch represents circular trajectories of the aircraft, where the radius of the circle decreases as the steering angle is increased. The saddle–node bifurcation at *δ*≈28° marks the loss of stability of the turn. In fact, as indicated by (*b*), for a steering angle *δ* past this point, the aircraft loses control without warning. Analysis of the state variables reveals that the tyres of the inner main gear saturate and the load is shifted to the outer main gear. The tyres of the outer main gear then saturate with an accompanying load redistribution to the nose and inner gears. The trajectory of the steering angle versus the forward velocity, past the saddle–node bifurcation point, also determines the final behaviour. For sufficiently large forward velocity, the aircraft may settle into a forward equilibrium (after a momentary LOC), but otherwise it may actually end up performing more complex motion, as depicted in figure 7*b*.

The stable branch (*c*) also represents circular trajectories of the aircraft but, compared with solutions on branch (*a*), the radius is much smaller. In other words, for the same steering angle and thrust settings, there are two different simultaneously stable turns with a hysteresis loop between them (as indicated by the vertical arrows). When the steering angle is decreased past the saddle–node bifurcation point at *δ*≈15°, the aircraft performs an outward spiral, as is shown by (*d*), until an equilibrium at a larger circular trajectory is found. In fact, the inner or the nose-gear tyres saturate, and then the particular tyres stay close to saturation, without any significant load shift to any of the other gears.

In the regions between the two Hopf bifurcation points, the equilibrium solution is unstable and stable periodic solutions and even chaotic motion can be found. Increasing the steering angle *δ* past the left-most Hopf bifurcation or decreasing it past the right-most Hopf bifurcation results in the onset of (initially small) stable oscillations. The aircraft will start shaking, giving the pilot a warning that the steering angle or thrust needs to be altered to move out of this region (figure 7*e*). When *δ* is brought into the region around 30°, then the periodic orbit is actually unstable and more complicated aircraft motion is found, which manifests itself as a complicated ground track (figure 7*b*).

The stable branch of equilibria in the lower right of the bifurcation diagram represents circular trajectories of the aircraft with a very small radius. The aircraft actually becomes stationary at very high steering angles, and (*f*) is an example. In this stopped position, a force impulse (e.g. a thrust impulse or wind gust) causes the nose gear to start dragging across the ground if no braking is applied (the aircraft will in fact stay stationary if the brakes are applied). The aircraft accelerates in a straight line, as is shown by (*g*), until a velocity on the top branch (*h*) is reached. This stable branch (*h*) of equilibria is somewhat counterintuitive in the sense that the radius of the circular trajectory decreases as the steering angle is decreased. Note that the radius of (*g*) is almost the same as the radius of (*a*). When the stability of this branch is lost by decreasing the steering angle *δ* to below the saddle–node bifurcation at *δ*≈50° the aircraft performs an inward spiral, as is shown by (*i*), until an equilibrium corresponding to a smaller turning circle is reached. Here, the inner gear tyres saturate, which is accompanied by a pitch-down motion. Then, the nose-gear tyres saturate, with a load shift to the outer gear.

A complete map can be constructed by conducting similar continuation runs at different thrust levels, as well as conducting two-parameter continuations of the saddle–node and Hopf bifurcations. This map can then be used to identify critical areas for further investigation, similar to that done at points (*b*) and (*d*). The reader can find further information on this procedure in Coetzee [38]. Additional studies have been conducted of a very large aircraft with a statically indeterminate gear [40], and of safe ground operating regions with the accompanying modes that lead to a LOC [41].

## 5. Outlook and conclusion

While use of bifurcation analysis is quite common in the analysis of the closed-loop response of military aircraft, it is less widespread in design and development of commercial aircraft. However, as was reported here, there have been some recent developments in this area.

Automation of flying vehicles is becoming more common, in particular in unmanned aerial vehicles (UAVs). The handling of highly flexible UAVs, such as X-HALE [42], will require sophisticated control algorithms; the scaling-up of platforms such as quad-copters will also result in highly nonlinear behaviour due to rotor dynamics. Bifurcation analysis is potentially of considerable benefit for the evaluation of these new aerospace applications, which are currently still quite poorly understood from the dynamical point of view. In particular, the insights gained by this approach will help enable the exploitation of nonlinearity, rather than avoid it as is traditionally the case. The same argument can also be made for more traditional vehicles such as future commercial airliners. Jones [43] reported the use of bifurcation analysis in a review of control challenges for future operations of commercial transport aircraft. Applications include using bifurcation diagrams of the form shown in figure 4 but with the addition of a fly-by-wire controller [44,45] to understand the requirements of automatic Return to Envelope Function strategies for aircraft recovery from upset. In addition, Kwatny *et al.* [46] used bifurcation analysis to conduct a rigorous assessment into the ability to control and to regulate the GTM in the vicinity of bifurcations, leading to the definition of safe sets in which the aircraft can be positively controlled.

The testing of fixed-gain and gain-scheduled controllers applied to the GTM has been conducted in Gill *et al.* [47]. To understand the sensitivity of the controllers, gain parameters were defined that scale the optimal gain values. It was argued that continuing in these gain parameters gives an understanding of the robustness of the controller and also enables the identification of the dominant control parameters. This might be thought of as a move towards understanding the effect of parameter uncertainty—a major challenge for the industry. Recent developments in this area relating to bifurcation analysis include the development of a continuation technique for the development of trajectories of stochastic differential equations [48] and the use of branch and bound methods [49]. The latter yields outer bounding sets of both the equilibrium and local bifurcation manifolds, and it was demonstrated with a study of the dynamics of a UAV with an uncertain centre of gravity location.

Industrial models are typically highly complex and may well contain hundreds to thousands of states and may include non-smooth characteristics captured via look-up tables or control saturations. This is in contrast to the models that are generally studied with bifurcation analysis, which might contain, say, from tens to a hundred states. To address this mismatch in model size, one of two approaches can be followed: either the large model is reduced without significant loss in terms of its faithfulness of description, or one needs to apply bifurcation analysis directly to the large-scale dynamical system. The latter approach is enabled by tools such as LOCA [50], but it has its own challenges and the learning curve is considerable. In addition, if one considers validation and certification, then a further challenge is to guarantee that all the possible nonlinear behaviours within the relevant parameter space are captured. A significant step in this direction is reported by Kolesnikov & Goman [51], who considered an industrial-scale model with non-smooth aerodynamics. With a matrix of output states, they solved the inverse dynamics to find solutions on all possible equilibria branches, which were then studied in detail by continuation.

From an industrial perspective, the challenges in introducing nonlinear dynamics into the engineer’s normal toolsets require demonstrating that these techniques provide added value over existing methods. Moreover, what is needed is a level of training and provision of user-friendly bifurcation software that seamlessly meshes to commercial software packages such as Siemens VLM or Matlab Simulink. Training aspects must cover both the formulation of the system so that it is suitable for bifurcation analyses, as well as the interpretation of the results—indeed, industrially relevant case studies are meant to assist in this. A level of intuition similar to that concerning, say, Bode diagrams needs to be developed by the practitioner for the interpretation of bifurcation diagrams. While many software tools are available, they were developed primarily for research purposes. The development of the DST is a first step towards the integration of bifurcation software into industrial tools. It is already helping companies such as Airbus solve industrially relevant problems.

## Authors' contributions

All authors collaborated on the reported work over many years; S.S. and E.B.C. provided models and industrial expertise; M.H.L., S.A.N. and B.K. led the bifurcation analysis and jointly supervised the different projects. All authors participated in writing the manuscript and gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

In addition to the financial support provided by Airbus, S.A.N. is currently supported by an EPSRC fellowship (EP/K0053737/1). Further support was provided by EPSRC CASE studentships, and support has been provided through the EPSRC programme grant ‘Engineering Nonlinearity’ (EP/K003836/1).

## Acknowledgements

The authors acknowledge the contributions from S. J. Gill and N. Terkovics, who kindly provided continuation plots for this paper, and from NASA, for the provision and support of the GTM.

## Footnotes

One contribution of 11 to a theme issue ‘A field guide to nonlinearity in structural dynamics’.

- Accepted February 26, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.