## Abstract

This article presents a systematic assessment of the use of numerical continuation and bifurcation techniques in investigating the nonlinear periodic behaviour of a teetering rotor operating in forward autorotation. The aim is to illustrate the potential of these tools in revealing complex blade dynamics, when used in combination (not necessarily at the same time) with physical testing. We show a simple procedure to promote understanding of an existing but not fully understood engineering instability problem, when uncertainties in the numerical modelling and constraints in the experimental testing are present. It is proposed that continuation and bifurcation methods can play a significant role in developing numerical and experimental techniques for studying the nonlinear dynamics not only for rotating blades but also for other engineering systems with uncertainties and constraints.

## 1. Introduction

The development of many accurate nonlinear numerical models of dynamical systems in engineering involves a level of experimental updating, ‘correcting’ or tuning, for example: the incorporation of parametric equations or experimentally derived correction factors. The extent to which these models are valid depends on the level of understanding of the behaviour of the physical system, which is—in turn—manifested in the modelling strategy and complexity together with any assumptions and experimental corrections.

In the aerospace sector, the use of nonlinear methods such as continuation and bifurcation tools is becoming more widespread. In particular, they are increasingly adopted to investigate nonlinear aircraft flight dynamics and control problems. However, the application of continuation and bifurcation methods has been limited to a small number of helicopter dynamical problems, such as flight mechanics [1–7], ground resonance [8,9] and examination of the rotor vortex ring state [10]. Furthermore, almost all of the investigations which use these nonlinear tools can be regarded as research studies, and it is still hard to find these tools widely adopted in industry for production aircraft. In recent years, the nonlinear aeroelastic stability of helicopter rotor blades was investigated by Rezgui *et al.* [11,12], using numerical continuation and bifurcation techniques. This investigation showed that these techniques are powerful in the identification of instability scenarios of rotor blades and uncovering the multiple solution structure driven by the nonlinearities in the rotor system. However, this work focused mainly on the applicability of the methods to investigate the effects of nonlinearities for standard helicopter blade stability problems. Nevertheless, this work led to the first practical application of the continuation and bifurcation methods for certification of production aircraft, which contributed towards the latest release-to-service of the AW159/Wildcat helicopter [13].

Unlike fixed wing aircraft, rotorcraft produce lift by spinning a number of blades about a rotor shaft, where each blade can constitute a separate dynamical system with its own degrees of freedom and hence its own behaviour. Creating an accurate numerical model for rotor dynamical analysis is very demanding due to the complexity of the aerodynamics and the structural properties of the rotor. For example, if the aerodynamics of the flow field is to be accurately represented, then all of the following characteristics need to be properly modelled: unsteady flows, compressibility, non-uniform inflow, reverse flow region, wakes, etc. Of course, this would substantially increase the model dimension. To carry out rapid blade stability analysis, however, it is customary to make use of a number of assumptions in order to reduce the size of the model to a practical level. This carries the risk of obtaining erroneous results, in particular in areas of fluid–structure interaction that have not been studied before, such as highly nonlinear regimes or novel rotor configurations.

On the other hand, the numerical tools used for bifurcation analysis can vary from simple brute-force time-integration techniques to the so-called numerical continuation methods. However, the accuracy and the correctness of the results produced by these tools are greatly dependent on the fidelity and validity of the nonlinear dynamical models used. In 2008, Sieber & Krauskopf [14] presented a novel experimental continuation technique that does not require an accurate numerical model. Instead, the physical model of the system under investigation is interrogated by the technique to produce the bifurcation diagrams, which describe its nonlinear behaviour. This experimental technique is often known as *control-based continuation* and can follow both stable and unstable solution branches (equilibria and periodic) as well as detect bifurcation points. There has been a number of studies in recent years that used the control-based continuation method for various purposes [15–20]. However, although these studies have already demonstrated how powerful the control-based continuation approach can be in tracing solution branches, without the need to know the underlying mathematical equations of the systems studied, the implementation of this method to more complex systems, such as fluid–structure interaction problems, can be challenging. In particular, if it difficult to achieve sufficient control of the required parameters because of constraints or limitations at the level of the physical system or the experimental set-up.

In recent years, the stability of rotor blades in autorotation was investigated using numerical continuation and bifurcation methods, and wind-tunnel testing. In [21], the coexistence of stable and unstable branches in the behaviour of autogyro rotors was demonstrated, both numerically and experimentally (not at the same time, as in the control-based continuation approach). The implications of such nonlinear behaviour on stability of autogyros from an engineering and flight safety perspective were discussed in [22]. In this paper, we revisit the not fully understood problem of autogyro blade flap–rotation instability, using experimental testing in conjunction with a numerical continuation and bifurcation analysis. The purpose is to highlight the benefits of the combined numerical–experimental approach in studying the dynamics of nonlinear fluid–structure problems. First, the blade flapping instability problem in autorotation is introduced with an overview of the principal challenges. Then, a description of the experimental apparatus and procedure are presented followed by a discussion of the obtained results. Finally, evaluation of the extent to which the continuation and bifurcation methods have helped in uncovering the physics underlying unstable blade behaviour is illustrated through a combination of numerical and experimental results.

## 2. Autogyros and instability of rotor blade in autorotation

Autorotation is a phenomenon whereby the rotor can sustain rotation relying only on the aerodynamic forces of the airflow passing through it. This process is exploited in aircraft, known as autogyros, to produce lift. For simplicity, the basic aerodynamics in autorotation can be viewed by considering the flow environment around one blade: in the general case, a certain section of the blade is absorbing energy from the airflow and the rest of the blade is adding energy to the flow. Hence, at the equilibrium condition, the net torque on the rotor shaft is zero and this constitutes steady autorotation. The autogyro achieves an autorotation state by slightly tilting the rotor disc backward allowing the air to flow upward and through the rotor blades as the aircraft moves forward (figure 1). To thrust the aircraft forward, an engine-driven propeller is typically used either in a tractor (puller) or pusher configuration [23].

Although autogyros were proved to provide significant advantages relative to other aircraft types, the greater operational envelope of the helicopter has relegated their use almost entirely to sport enthusiasts; hence they are less familiar to the wider public. The autogyro fraternity knows, through experience, that a rotor operating in high-speed edgewise flight can exhibit rather undesirable characteristics in the form of blade flapping instability or rapid rotor speed decay, especially if it is lightly loaded (low thrust coefficient). However, the mechanism governing this potentially very dangerous behaviour is poorly understood and a better understanding of the lifting rotor dynamics that can lead to such wayward behaviour is required.

The lack of understanding of autogyro stability together with inadequately trained pilots resulted in a number of fatal accidents. In an attempt to categorize the autogyro crashes on the basis of their causes, the Popular Rotorcraft Association identified aircraft instability—where the aircraft becomes unstable and uncontrollable—as one of the two major contributions to the accident record [24]. Violent blade flapping leading to in-flight rotor blades striking other aircraft components is one manifestation of this, along with the loss of rotor speed. These outcomes are directly related to the rotor blade stability and are the focus of this paper.

### (a) Blade flapping in autorotation

Like in helicopters, the blades of an autogyro rotor are allowed to flap to avoid the problems of lift asymmetry. The flapping behaviour of the blade in forward flight will create an oscillatory motion, in a way that the peak flapping amplitude of the oscillation is achieved over the nose of the aircraft and the minimum is achieved over the rear of the aircraft. This results in the rotor tip-path plane being tilted back as viewed from the side of the aircraft. This rotor disc tilt increases with forward speed [23], which in turn increases the flow rate through the rotor. This means that rotor shaft angle (the hub plane angle) can be reduced even further to keep steady-level flight. In forward flight, the aerodynamic forces provide excitation to the flapping blade, primarily at once per revolution and hence establish a periodic forcing in the rotor system.

The flapping oscillatory motion of the blade also induces cyclic forces in the plane of rotation, owing to the conservation of angular momentum of the blade. These forces are called Coriolis forces and they require the incorporation of lead-lag hinges at the blade root to allow blade motion in the plane of rotation to eliminate structural fatigue. Finally, as well as the aerodynamic and Coriolis forces, there are the centrifugal and gravitational forces acting on the blade.

### (b) Stability of autorotation state

It was seen that autorotation is a phenomenon that relies on balancing the aerodynamic torque acting along the blade. Hence, it is logical to investigate whether an autorotation state is always achievable. Of course, if it cannot be sustained during flight, then this will be regarded as a form of instability. The stability of the autorotation state was first investigated in axial descent by Wimperis [25] in 1920. He showed that the autorotational state of a rotating blade section is dependent on the local inflow angle and the local angle of attack, as well as the section’s aerofoil characteristics manifested in the lift and drag coefficients. This was investigated using the so-called *autorotational diagram* and it was shown that both the accelerating and decelerating forces acting on the blade section are stabilizing. It was also shown that there is a maximum pitch angle above which autorotation is not possible, regardless of the inflow angle value. This represents the stall condition, which causes only a decelerating force to exist. However, although the autorotational diagram can be used for a single blade section, constructing a similar diagram for the whole rotor blade is very complicated. This is because autorotation equilibrium in this case is determined by the cumulative effects of the forces and flow velocities acting along the blade. Similarly, constructing an autorotation diagram for the forward flight case is also challenging due to the oscillatory behaviour of the blade and indeed for any blade section.

Nikolsky & Seckel [26] extended the above analysis to rotors in axial autorotation. They considered the effect of stalling on the stability of autorotation. They showed that when stall effects are included, two trim solutions can be found. The first represents the normal stable autorotation state, while the other is an unstable autorotation condition. The analysis also revealed that for small blade incidences, the stability of the blade is evident. However, for higher blade incidences, there is a risk that flow disturbances can cause the blades to stall, because of the unstable trim points being close to the normal autorotation state. Nikolsky and Seckel also illustrated that there is a maximum incidence angle (about 8.8° for the example rotor used) above which axial autorotation cannot exist.

### (c) Blade stability in autorotation

As far back as 1936, Wheatley [27] found that autogyro rotor lead-lag oscillations are mainly the direct effects of blade flapping (Coriolis forces) and that the influence of the aerodynamic forces in the plane of rotation are secondary. Understanding the lead-lag motion of an autorotating rotor can help in knowing the amount of damping required by the blade. The analysis was also endorsed by flight test validations. In 1978, Wei & Peters [28] used perturbation and Floquet methods to study the blade lag stability in autorotation, including the effects of advance ratio (ratio of aircraft forward speed to blade tip speed), blade elastic coupling, blade in-plane natural frequency and rotor trim condition. This analysis revealed that blades operating in the autorotation condition are considerably less stable compared to the powered flight regime. This study was, however, limited to cases of advance ratios less than 0.5.

Floros & Johnson [29,30] investigated aspects of blade flap-lag stability and found that it was difficult to find trim solutions above an advance ratio of 2. They also identified that the trim solution in autorotation is not unique, which raises the question whether a manoeuvre could cause the rotor to change abruptly between different states. Rigsby [31] also used Floquet analysis to investigate the stability of an autorotating rotor. Although the results identified some areas where the stability is reduced, the analysis did not predict any form of blade instability.

Another source of rotor instability in autogyros is the interaction between the blade elastic forces with the aerodynamic and inertial forces. Although the aeroelastic effects are not the subject of this paper, it is worth mentioning that there is a strong coupling between the rotor speed and the blade degrees of freedom in bending and twist [32,33], and that blade elasticity can drastically reduce the rotor stability [29].

Finally, there is extremely limited experimental research in the available literature that addresses blade stability for autogyro rotors, most testing focusing instead on performance and flying qualities. De Silva [34] performed wind tunnel tests to quantify the extent of flapping at high advance ratios. The tests showed that as the advance ratio increases, the amplitude of the blade flapping angles increases first linearly and then exponentially, depending on the rotor disc loading. However, the results of the experiments did not distinguish between real flapping instability and high flapping angles. On the other hand, Wheatley & Bioletti [35] confirmed one of the basic concepts about autorotation stability by conducting wind tunnel tests using a 10 foot diameter rotor. The tests showed that there is a critical angle of attack (rotor tilt) below which the rotor can self autorotate and above which the rotor would rotate in the reverse direction. However, since the rotor had only one degree of freedom, azimuthal rotation, then it was inevitable to spin either in autorotation or in a reverse direction.

## 3. Continuation and bifurcation methods for rotating blades

The basic idea of numerical continuation and bifurcation techniques is the calculation of the steady solutions of a dynamical system as one of its parameters, called the continuation parameter, is varied across a pre-defined range. The computed solutions construct a number of branches that could be either stable or unstable. To determine the stability, either an eigen or Floquet analysis is carried out at each computed solution, depending on the nature of the solution. For instance, in hover the blade behaviour can be considered to be in equilibrium (fixed points), hence an eigen analysis is carried out for stability, whereas in forward flight, the blades behave in a periodic manner (limit cycles) due to the rotor lift asymmetry, and hence Floquet theory is used to determine the stability.

A bifurcation is the qualitative change in the system behaviour as a parameter is varied. In other words, when the stability of a system is changed or lost, the system bifurcates. The points at which these stability changes happen are called bifurcation points. When the system is nonlinear, new solution branches may emerge from the bifurcation points, leading to the presence of multiple solutions for the same set of system parameters. The identification of these different solution branches helps to uncover the global dynamics of the system. Of particular interest is when the blades, for example, are locally stable for small disturbances but not necessarily for large ones, and vice versa.

Therefore, the strategy in implementing continuation and bifurcation methods is to follow one solution branch as one or more parameters are varied to locate bifurcation points. The emerging branches are then followed to construct a more complete picture of the system dynamics (bifurcation diagram). Furthermore, other advantages of continuation methods, compared with other time history or frequency domain methods, are their efficiency and accuracy in following the solution branches as well as in detecting and identifying the bifurcation points. The different types of bifurcations that can occur in equilibria or periodic orbits are not discussed in this paper; the reader is referred to general texts such as [36] for more background on the subject.

There are several freely available continuation software with different levels of maturity and robustness. AUTO [37], Continuation Core (COCO) [38], MATCONT [39] and CL_MATCONT [39] are examples of the most widely used packages. In this analysis, the continuation and bifurcation software AUTO was used. AUTO is open source software for continuation and bifurcation problems of ordinary differential equations, originally developed by Eusebius Doedel, with subsequent major contributions by several people and it is currently available on a number of platforms [37,40]. Besides many other types of equations, AUTO can perform extensive bifurcation analysis of ordinary differential equations of the form
3.1subject to initial conditions, boundary conditions and integral constraints. Here *x* is the state vector and *p* denotes one or more parameters. *n* and *m* are the numbers of states and parameters, respectively. Equation (3.1) is written in the generic (nonlinear) state-space form, where the state-derivatives are functions of the states and some parameters.

The main type of steady solution which describes the rotor blade behaviour in the conventional forward flight operating envelope is a periodic orbit (limit cycle solutions). Conventionally, helicopter rotor models have a constant rotor speed and are written in the non-autonomous form, in that the independent variable *t* appears explicitly in the equations. In fact, the blade azimuth angle *ψ* is a non-dimensional form of the time variable *t*. Unlike in helicopters, the rotor in forward autorotation has a variable rotational speed. Hence the rotor is not forced to rotate at a fixed frequency and the blade azimuthal angle *ψ* would then need to be modelled as a state variable. The rotor in forward autorotation is therefore a self-excited dynamical system.

## 4. Description of experiment

The experiments were performed in the University of Bristol low-speed open jet wind tunnel (figure 2*a*). The wind tunnel is a closed return system with a 1.68 m long open working section. The diameter of the jet is 3 ft 6 in (1.1 m) and the maximum attainable velocity is about 33 m s^{−1}. The experimental rig (figure 2*b*) comprises a two bladed teeter rotor of 1 m in diameter. The blades are rigidly connected and free to flap about a hinge located at the shaft axis. The skeleton of the rotor rig, including the hub system, is modified from a radio controlled helicopter. The flapping angle of the blades is measured using a magnetic encoder which has a resolution of 0.15°±0.07°. The pitch of each blade is monitored using different magnetic encoders that are connected to each blade via a pulley and tooth belt system, resulting in a total pitch resolution also of 0.15°±0.07°. The signals from these three sensors are transmitted via ZigBee wireless telemetry to a personal computer.

The rotor speed and azimuthal position are monitored using an optical encoder connected to the rotor shaft. This encoder is connected to the PC via a dSPACE interface (dSPACE Ltd, http://www.dspace.ltd.uk, accessed February 2015). A six component load cell (JR3-100M40A (JR3 Inc., http://jr3.com/info/about.php, accessed February 2015)) fitted below the fuselage is used to measure the forces and moments acting on the rotor rig. The airframe is designed to represent a scaled version of the Magni VPM-16 autogyro, modified by having a closed cockpit. The purpose of this airframe is to cover the components of the rotor support frame and provide a smoother aerodynamic shape, and hence was not designed to match the scale of the rotor dimensions. Finally, a safety mechanism is implemented to prevent damage to the rotor rig should the blade motions become unstable. This mechanism is located in the rotor hub and uses two spring loaded valves, which are released by a mechanical locking system when contacted by a blade exceeding the maximum allowed flapping angle of 23°.

The procedure was to monitor the rotor speed and blade flapping angles in autorotation, as the wind speed (*U*), shaft angle (*θ*_{shaft}) or blade collective pitch angle (*θ*_{col}) were varied. The first series of tests were conducted at a collective pitch angle of 1°. In each test, the shaft angle was inclined at a chosen value, where the blade flapping angle and rotational speed were recorded as the rotor was subjected to different wind speeds, taken at 1 m s^{−1} intervals. The aerodynamic forces and moments were also measured. At every wind speed value, the rotor was allowed enough time to settle down to steady autorotation, then the readings for the flapping angle and rotational speed were recorded for 30 s and the readings for the aerodynamic forces were recorded for 6 s. These experiments were repeated twice to ensure reproducibility of the results. The rotor often required to be pre-rotated (by hand or using string wound around the shaft) for autorotation to take place. This set of experiments was repeated for collective pitch angles of −1°, 0°, 2° and 3°.

It is worth noting that some measurements were not possible due to unforeseen rotor shaft vibration (shaft resonance), which often occurred when the rotor rotated close to average speeds of 500 r.p.m. (±1%), 960 r.p.m. (±1%) and 1270 r.p.m. (±1%). The last two frequencies were the most potentially damaging ones for the rotor and hence it was not practical to keep the rotor spinning at those speeds if excessive vibration occurred. Operating the rig close to these frequencies also affected the accuracy of the measurement, which resulted in the largest discrepancies in rotor speed and flap angle. Excluding those cases, the maximum scatter in the measured rotor speed was about 35 r.p.m. Whereas, the maximum scatter in the flapping angle amplitudes was about 1.5°.

The next step of the experiment was to identify, within the wind speed range, any areas where the blade dynamics change. This was done in two stages. Firstly, by investigating the rotor behaviour at the high and low limits of achievable wind speeds. Secondly, by perturbing the rotor velocity during the normal autorotation regimes. This was simply done by introducing a level of friction for a short period of time to slow the rotor down.

## 5. Experimental bifurcation analysis

Figure 3 depicts the variation of the rotational speed and flapping angle at steady autorotation state. Since the recorded flapping data were periodic with rotor azimuth position, the average peak values for each cycle were computed over the period of the data recording. These peak values of the oscillating flapping angle are plotted in figure 3 for different tunnel speed values and rotor shaft angles. For the rotor velocity, the averaged mean values were used instead, in order to minimize the effects of the once per revolution contributions. In steady autorotation, the amplitude of the rotor velocity oscillation (lead-lag rate) is very small, so peak and mean values can be used interchangeably without affecting the analysis, particularly when comparing with the numerical results. The best-fit curves for the measured data were also plotted. These curves were extrapolated to obtain an illustration of the rotor speed and peak flapping angle variations at higher wind speeds. Figure 3 illustrates that the rotational velocity increased almost linearly with the forward wind speed, while the flapping angle had an inverse relationship with the wind speed. This shows that the more flow that went through the rotor the faster it rotated but at the same time it affected the flow environment in a way that forced the steady-state flapping angle amplitude to reduce for higher speeds.

The experiment also revealed that at each setting of shaft angle, autorotation was not possible below a certain wind speed value, which indicated the presence of an instability point. In other words, the presence of a bifurcation point. Instead, regardless of how high (within the permissible range) or low the initial value of the rotor speed was, the blades’ rotational speed always decayed and the flapping oscillation diverged until the blades made contact with the safety mechanism, which was then activated (figure 4). It should be noted that it was difficult to precisely determine the minimum autorotation wind speed, particularly with the adverse effects of tunnel flow perturbations: the instability points illustrated in figure 3 were located within a ±1.5 m s^{−1} band of the stated wind speed value.

The subsequent experiment was to test the ability of the rotor to autorotate at flow speeds higher than the minimum autorotation airspeed, but in this case starting from different low initial rotor velocities at each flow speed. To achieve a level of control in the initial rotor speed values, the rotor is first operated in stable autorotation after setting the appropriate wind tunnel speed. Next, a level of frictional torque is applied at the rotor shaft to reduce the rotor speed down to the desired value. The top of the wireless telemetry shroud was found to be the best place to apply this resisting torque. For every wind speed and shaft angle setting, a rotor speed threshold was found, above which the rotor will be able to achieve steady autorotation. Figure 4 depicts the results for different attempts of initial rotor speed conditions when the shaft angle was set to 7°. By repeating the tests a number of times, it was possible to define a rotor speed boundary that separates the two rotor scenarios (presented by a red dashed line in figure 4). At this boundary, the behaviour of the rotor appears to be steady but due to flow disturbances it will move away from this point either to the stable autorotation state or the blade flapping diverges and activates the safety mechanism. This rotor speed boundary together with the corresponding peak flapping angle variation represent an unstable periodic branch of a rotor dynamics in autorotation.

The diagrams in figure 4 resemble experimental bifurcation diagrams of the physical teetering rotor. They show the presence of stable and unstable autorotation branches, which get connected at the bifurcation point. Because of the latter and also because of the unidirectional transient behaviour of the blades at speeds below the bifurcation point wind speed, the instability point resembles a *fold bifurcation point*: it defines the minimum wind speed value for which steady autorotation will exist. Moreover, the approach used to find the above bifurcation diagrams is typically known as a *brute-force bifurcation* method.

Furthermore, from the above experiments, two other issues need to be clarified. The first is related to the behaviour of the rotor at a wind speed below that of the bifurcation point. It seemed during the experiments that all attempts to operate the rotor in these low wind speed cases ended up by activating the safety mechanism due to the large flapping angle before the rotor speed decayed to a halt and not the opposite way around. This could mean that the rotor at low speed might be attracted to another stable periodic branch that resides beyond the limit of the maximum flapping angle. Of course, if that is the case then the system features important multi-attractor characteristics and hence continuation and bifurcation methods are really an essential tool for investigating its dynamics, when the appropriate numerical rotor model is available.

The second issue concerns the location of the bifurcation point as the shaft angle was changed. Figure 5 is a two-parameter experimental bifurcation diagram and shows that the wind speed value at which the fold bifurcation occurs increases as the shaft angle is reduced. This implies that the fold bifurcation points also specify the minimum permissible shaft angle required for autorotation at those wind speed values. In order to clarify the effects of the above results on autogyro flight performance, the average aerodynamic lift values generated by the rotor are superimposed as contours in the same figure. Figure 5 illustrates that the fold bifurcation line lies in an area where the rotor produces low values of lift (below 1 N to just above 3 N). This shows that when the rotor is unloaded by reducing the shaft angle or the flow wind speed, its stability will be determined by whether or not the fold bifurcation line is crossed. Furthermore, the higher the wind speed the higher the minimum lift value generated by the rotor. This means that the minimum permissible lift value increases at higher wind speed. This is an important result because, for high-speed flight, the autogyro pilot has to reduce the shaft angle to maintain the lift and weight balance or to further off-load the rotor for specific manoeuvres. Thus, if the shaft angle is reduced enough so that the fold bifurcation point is reached then the rotor will not be able to sustain stable autorotation and, hence, the blades will start to flap in a oscillatory divergent manner.

Finally, the effects of varying the collective pitch angle were investigated. The shaft angle chosen for this case was *θ*_{shaft}=7°. Attempts to get the rotor to autorotate at *θ*_{col}≤−2° and *θ*_{col}≥4° did not succeed over the available speed range of the tunnel. Figure 6 depicts the variation of rotor lift in the blade pitch angle and the wind speed parameter plane. In this case, the fold bifurcation line still lies in an area where the rotor produces low values of lift but the range of values is relatively bigger (below 0 N to just above 7 N) than the range shown in figure 5. Furthermore, it can be seen that in the range of 11 and 17.5 m s^{−1}, two bifurcation points at every wind speed exist; one at a small collective pitch angle and another at a higher angle. For example, at *U*=16 m s^{−1}, two fold bifurcation points exist at approximately *θ*_{col}=−0.9° and 2.5°. From this plot, one can easily visualize the appropriate bifurcation diagram for the case of *U*=16 m s^{−1}, when the pitch angle is taken as the continuation parameter. Figure 7 depicts this bifurcation diagram. For *θ*_{col}=0° and 2°, the mean rotor speeds and the peak flapping angles at the unstable branch were found in a similar way to that used in figure 4. Figure 7 clearly illustrates that autorotation is only possible within a range of blade pitch angle, which is dependent on wind speed and shaft angle. At the limits of this range, there are two bifurcations that define the points after which the torque balance in the rotor can no longer be achieved. This type of bifurcation diagram is usually called an *isola*. The above results represent a new qualitative interpretation of autorotating rotor stability. The multi-attractor nature of nonlinear dynamical systems suggests that the instabilities identified above could lead to new attractors that will constitute a more global view of the system behaviour than has previously been available.

## 6. Description of rotor model

It is apparent that the rotor behaviour is very nonlinear and that analysing such a system numerically necessitates the use of nonlinear analysis methods with a numerical model that can capture the fold bifurcation seen in the wind tunnel tests. What makes this type of investigation even more difficult is the fact that the number of parameters involved together with the different flight manoeuvres and rotor types that need to be considered is large. From an engineering perspective, it is important to be able to perform numerous quick analyses in order to minimize the problem dimension to a manageable level. For this reason, it was necessary in this investigation to create a low-order nonlinear rotor model that nevertheless retains acceptable fidelity in the aerodynamic and dynamic representations. Since the focus here is on flap–rotation instability (coupled flapping motion with blade rotation about the main rotor shaft)—as opposed to, say, blade divergence and flutter or blade/pitch-link coupling dynamics—it is adequate to assume rigid blades and to consider degrees of freedom both in the plane of rotation and normal to it. The omission of blade elastic characteristics reduces complexity and allows easier insight into the blade dynamics.

The model developed for this study represents the dynamics of a two-bladed rigid teetering rotor. It is in the nonlinear continuous state-space form where is the state vector, is a vector of *m* parameters (e.g. blade pitch angle, wind speed, etc.) and *f* is the set of *n* nonlinear equations (equations of motion, inflow model, etc.). The formulation as a set of nonlinear continuous ordinary differential equations makes it amenable to the application of bifurcation analysis. The following simplifying assumptions—demonstrated to provide reasonable model performance output for autorotating rotors [41]—are adopted in the aerodynamic representation:

—

*Two-dimensional steady aerodynamics*. A two-dimensional individual blade element approach is adopted. The lift, drag and pitching moment of each element are calculated numerically using nonlinear look-up tables as functions of angle of attack and Reynolds number. Experimental data for a NACA0015 aerofoil section are used for each blade element [42]; since it is of interest to investigate the high advance ratio condition, in which reverse flow is expected to occur over a portion of the rotor disc, the look-up tables cover the whole 360° range of angle of attack.—

*Unsteady flow and vortex shedding effects ignored*. While unsteady flow characteristics such as dynamic stall clearly play an important role in rotor aerodynamics, the instability mechanisms of interest here are expected to exist even when these are neglected. This assumption was made in the interest of reducing model complexity and hence computation times.—

*Compressibility effects neglected*. Because the blade tip speed is generally low in autorotation even at high forward speed due to low rotor speed, compressibility effects were neglected but they can be easily incorporated in the look-up tables.—

*Blade/blade and blade/hub interactions ignored*. Again, these were neglected for the purpose of simplicity and efficiency in generating solutions.—

*Tip losses approximated*. It is assumed that the lift forces at the tip and blade root reduce to zero. For this analysis, a tip loss factor of 97% is used; other tip losses are ignored.

Rather than describing flapping by its Fourier components in the non-rotating frame, the individual flapping coordinate, *β*, is retained for each blade and also its azimuthal position, *ψ*. For the teeter rotor, both blades are assumed to be rigidly connected and hence have one flapping degree of freedom; also, blade rotation together with blade lead-lag motion are represented by the same state variable, namely the rotational degree of freedom around the shaft axis, *ψ*. The flapping coordinate for the blades, *β*, is dependent on *ψ*. The differential equations of motion for the rotor in both the flapping and the rotation senses are second order, giving a total of four rotor states (*ψ*, , *β* and ), which are allocated for each rotor blade according to table 1. A compact form of the model equations is presented in appendix A, while the full derivation of the equations can be found in [43].

The induced flow through the rotor is captured via a three-state Pitt–Peters dynamic wake model [44–47]. This representation forces a linear distribution of the induced velocity along the rotor disc radius and a sinusoidal variation along its azimuth. The model was originally developed for powered helicopter rotors, and modified by Houston for rotors in autorotation [48–50].

Scrutiny of the overall model equations shows that the model contains geometric, inertial and aerodynamic nonlinearities and bifurcation analysis is an ideal tool for the study of the system stability and its dependence on parameter variations. The model is sufficiently smooth, despite the look-up tables, for implementation in numerical continuation software; AUTO has proved to be both versatile and efficient in this work. The model was implemented in Matlab, and AUTO was run from within this environment.

The model is seventh order, the state vector being . A number of the model parameters can be used as continuation parameters. However, to understand the rotor stability during operational conditions, the main parameters of concern are flow speed (*U*), blade collective pitch angle (*θ*_{col}) and rotor shaft tilt angle relative to the flow (*θ*_{shaft}); the latter is subsequently referred to as simply ‘shaft angle’.

## 7. Numerical analysis

The experimental results presented in §5 clearly demonstrate that the blade dynamics in autorotation is nonlinear. The obtained experimental bifurcation diagrams illustrate that the loss of instability for lightly loaded rotors can be understood as the consequence of operating near fold bifurcations. In order to understand these experimental results, particularly beyond the experiment limitations (for example, higher wind speed or higher flapping angles), a numerical analysis was carried out. The appropriate physical rotor properties were adapted in the numerical rotor model. The primary objective of this investigation is to evaluate the numerical model correlation to the experimental results, in particular its ability to capture the fold bifurcation points of the rotor.

### (a) Model tuning

The blade dimensions, mass and aerofoil aerodynamic characteristics of the experimental autogyro blade were used in the numerical model. Initial spot checks using time history simulations predicted rotor speed values much greater than those values achieved in wind tunnel tests. The main reason for this large difference in rotor speed is believed to be the simple aerodynamic representation used in the model. This discrepancy in rotor speed also led to dissimilarities in other aspects of the model, notably the rotor flapping angle and forces produced by the rotor. Further scrutiny of the results revealed that the high value of the predicted autorotational speed is related to the poor estimation of the rotor aerodynamic torque, which caused the torque balance to be realized at a higher rotor speed than that experimentally measured.

Therefore, it was necessary to correct for some of the aerodynamic characteristics that were not accounted for or were not accurately modelled. For instance, accurate blade tip and root losses, blade-to-blade interaction, unsteady aerodynamics including dynamic stall, airframe and tunnel interaction effects, inaccurate rotor downwash, etc. For simplicity, it was decided to use a friction term in the rotation sense of the rotor to correct for those mis-represented characteristics. The use of this crude assumption was not expected to produce an exact match between the numerical and the physical model for all the rotor states. However, the rationale was that this correction would narrow the discrepancy gap and produce a reasonable match in rotor speed, at the same operating parameters, without increasing the complexity of the model. Furthermore, in reality, an amount of resisting torque was present in the rotor shaft due to friction in the hub bearings and the swash plate.

The frictional torque correction term is assumed proportional to the rotational blade velocity:
7.1where *ζ*=*f*_{1}(*θ*_{shaft})+*f*_{2}(*θ*_{col}) is a friction ratio and is a function of both the shaft and the blade collective pitch angles. The following functions for *f*_{1} and *f*_{2} were found to provide a reasonable fit to the experimental data:
7.2and
7.3where the shaft and blade pitch angles in the above equations are in degrees. The tuning of the friction coefficient *ζ* was based on trimming the rotor model by only matching the peak rotor velocities to those found in the experiment at the stable autorotation branches. It should be noted that it was not possible to trim the numerical model at low values of forward speed (below 25 m s^{−1} for shaft angle of 7° and pitch of 1°). However, this was resolved by extrapolating the experimental rotor speed values above the 30 m s^{−1} limit of the wind tunnel. As will be discussed in the following section, the reason for not achieving trim at low flow speeds was because the location of the fold point in the numerical model was sensitive to the torque coefficient. Finally, the nonlinear least-square curve method in Matlab (lsqnonlin) was used to fit the torque coefficient data gathered to the function described in equations (7.2) and (7.3). Furthermore, one has to bear in mind that even if the rotor speed is well tuned, the other states are not expected to be closely matched because of this crude experimental correction. However, the aim was to be able to produce at least a qualitative agreement in the topology of the bifurcation diagram.

### (b) Bifurcation analysis

Figure 8 shows two projections of the bifurcation diagram generated using the tuned model, with wind speed as continuation parameter. Note that all solutions are limit cycles, for which the maximum amplitude is shown for each state variable component (rotor velocity and flapping angle). The continuation started from a periodic orbit at wind speed of 60 m s^{−1}, which is then reduced to trace out the first stable autorotation branch. As shown in figure 8, the end of this branch is a fold bifurcation point which is the birth of an unstable branch. By increasing the wind speed a second fold point is found, where in this case the branch turns back and becomes stable once again. This second autorotation branch is characterized by higher flapping angle amplitude and slower rotational velocity than the first stable autorotation region. As the wind speed is decreased this second autorotation branch loses stability at a *branch point bifurcation* and two new unstable periodic branches are formed. The flapping angle graph shows this latter result more clearly than the graph of rotor velocity. Because of the bifurcation diagram geometry near this point, it is usually called a *pitch-fork bifurcation*, and it is associated with the breakdown of symmetry in the dynamical system. At slightly increased wind speed, both branches bifurcate at new fold bifurcation points resulting in two new stable branches (branches 3 and 4; figure 8).

Like the experimental bifurcation diagram of figure 4, the above results confirm the presence of the fold bifurcation, although located at higher wind speed values. These results illustrate that the qualitative stability characteristics of a rotor in autorotation can be predicted, when continuation and bifurcation techniques are adapted, even though the numerical model is relatively simple and of low order. Furthermore, similar results (not shown here) confirm that the wind speed value at which the fold bifurcation points exist increases as the shaft angle is reduced. Therefore, the boundary of rotor stability can also be constructed by plotting the fold bifurcation line in the two parameter space of the wind speed and the shaft angle [22].

Figure 8 shows that in the wind speed range of 27–56 m s^{−1} the rotor can autorotate at two different rotational speeds. This means that if the rotor speed starts initially from a lower value than that of the unstable branch, then the rotor states will get attracted towards the second stable branch and the rotor should be able to achieve autorotation, but at lower rotational speed. However, the problem with the latter is the high flapping oscillation (at least 18° in amplitude). This could be the reason why in the experimental results in §5 there was no sign of such rotor behaviour since the safety mechanism will not allow it to happen. Furthermore, the existence of a second fold bifurcation point means that the rotor can undergo hysteresis-like behaviour. This happens when the rotor is initially autorotating along the second stable branch and if the wind speed is increased above 56 m s^{−1} the rotor velocity will jump from about 1300 r.p.m. to over 3000 r.p.m. (see right green arrow in figure 8).

In figure 8, the second stable branch bifurcates to effectively two new stable branches (the parameter range of unstable branches 3 and 4 is very small and it would be very difficult to achieve any experimental match). The new stable branches (3 and 4) extend from a wind speed value of 36 m s^{−1} to less than 5 m s^{−1}. Interestingly, even though these two branches trace different flapping angle amplitudes, they have the same rotational velocity. The flapping angle and the rotational speed graphs for blade 1 at branch 3 are exactly the same as those of blade 2 at branch 4 lagged by 180°. Similarly, the flapping angle and the rotational speed graphs for blade 2 at branch 3 are exactly the same as those of blade 1 at branch 4 lagged by 180°. What is more interesting is that, for the same branch, blade 2 does not repeat the same behaviour as blade 1 after 180° as is the case in a normal autorotation regime. This is proof of the loss of symmetry between the blades of the teetering rotor and as previously mentioned, it is a characteristic of the pitch-fork bifurcation point.

Finally, to determine whether the new tuned model can also provide similar behaviour to that of wind tunnel test when the blade collective pitch angle is varied, a continuation run was performed for *θ*_{shaft}=7° and *U*=30 m s^{−1}, taking the collective pitch as the continuation parameter. Figure 9 depicts that this results in a qualitatively similar bifurcation diagram to that obtained using the experimental measurements (figure 7). The results clearly confirm that autorotation is only possible within a range of blade pitch angle, the limits of which depend on wind speed and shaft angle. At the limits of this range, there are two bifurcations that define the points beyond which autorotation is impossible. This bifurcation diagram also confirms that the stable and unstable branches constitute an *isola*. Note that the stable autorotation branch pertains to branch 1 of figure 8.

It can be seen in figure 9 that the maximum collective pitch angle for autorotation to exist is less than 2°. This obviously is less than expected for this rotor and it is thought that the reason for this low value is the simple tuning procedure that was implemented in this analysis. In general, both experimental and numerical results confirmed the dependence of fold bifurcation on collective pitch. Blade twist as well as blade aerofoil shape are also parameters that can influence the blade stability and, hence, can be considered for design of more stable rotor blades in autorotation.

## 8. Extended experimental analysis

In order to experimentally confirm the presence of a second stable branch, a wind tunnel test was conducted. This time, the collective pitch angle was reduced to 0° since figure 7 suggests that the rotor would be able to rotate at slower speed and smaller flapping angle than the maximum allowed limit. Interestingly, when the shaft angle was set to 7° the rotor was able to autorotate at about 100 r.p.m. at a wind speed of 27 m s^{−1}. Moreover, as the wind tunnel speed was very slowly increased to about 28.5 m s^{−1}, the rotational speed increased slightly and in a steady manner to 120 r.p.m., which proves the presence of a second small stable branch. However, if the wind speed is raised to 29 m s^{−1}, the rotor can no longer sustain autorotation on this branch; instead, the rotational speed jumps very rapidly to over 950 r.p.m. These results are plotted in figure 10*a*. The wind tunnel test sweep conducted at *θ*_{col}=0° shown in figure 10*b* confirms that the branch characterized by the higher rotor speed value is the main stable autorotation branch. Note that in the gradual wind tunnel sweep, it was not possible to get the rotor to autorotate at flow speeds above 22 m s^{−1} because of early onset of excessive vibrations. However, the rotor managed to autorotate at 29 m s^{−1} when the transient behaviour started from a relatively very low initial rotor speed condition.

The flapping amplitudes at branches 3 and 4 in figure 8 are very high. This makes it difficult to experimentally find a similar rotor behaviour due to the physical flapping limitations incorporated in the rotor rig. However, with the collective pitch at 0° and the shaft angle set to 7°, the blade starts rotating in an antisymmetric manner at lower speeds than in figure 10. For example, in figure 11, where the wind speed is 25 m s^{−1}, the average of the flapping oscillation is no longer 0°: instead, it is about 7°. This means that the blades are not repeating the same flapping motion, as they do in normal autorotation (where each blade would be at the same flapping angle for the same azimuthal position, i.e. average flap amplitude zero).

## 9. Conclusion

An evaluation of autogyro rotor stability by numerical continuation and bifurcation tools and using wind tunnel experimentation has been undertaken. The case study of rotor blade instability in autorotation was used to discuss the real engineering issues in predicting unstable behaviour of complex nonlinear aero-mechanical systems. It was shown in this example that obtaining a valid rotor model is an essential part of the investigation, although the continuation technique provided a powerful tool to ease decision-making in arriving at an improved low-order model during the analysis.

The systematic procedure suggested by this work involves a collaborative use of continuation and bifurcation analysis, using a low-order but representative model, with experimental testing (not necessarily at the same time, as in the control-based continuation approach). Assuming that bifurcation diagrams for a corresponding physical system can be produced experimentally, efficient use can be made of the experimental process to supplement and verify the numerical results as well as to tune the numerical model. It was shown for the autogryo rotor system that prediction of attractors beyond the maximum flapping angle was possible and that the numerical results are in good agreement with those generated from experiments—at least qualitatively in the highly nonlinear region. The approach was able to identify new blade behaviour arising from the nonlinearity of the system, namely asymmetric autorotation branches following a pitch-fork bifurcation of the periodic orbits (also verified experimentally). These have not previously been characterized in this way.

The continuation method detected fold bifurcation and branch points, which cannot be achieved by other methods. Indeed, some approaches may mistakenly suggest unstable solutions in the region where there should not be any solution, such as beyond a fold point where there is no local solution branch. As demonstrated for the rotor system, the combined numerical–experimental procedure facilitates the construction of full bifurcation diagrams, even where state or parameter constraints exist. In this way, a comprehensive insight into the multi-attractor nature of the autorotating blade behaviour was obtained; this provided a new explanation of autogyro rotor instability. Although not shown here, a further advantage of the continuation-based approach is that information from several bifurcation diagrams can be merged to produce a three-dimensional bifurcation diagram. This can be used to establish stability boundaries in parameter space; alternatively, two-parameter continuation can generate the stability regions.

## Competing interests

We declare we have no competing interests.

## Funding

The authors would like to thank the Algerian Ministry for Higher Education and Scientific Research for funding the work presented in this paper.

## Acknowledgements

The authors are also very grateful to Mr Peter C. Bunniss from the University of Bristol for his continuous support and guidance throughout the undertaking of the work presented in this article.

## Appendix A. Rotor model equations

The full derivation of the equations of motion can be found in [43]. However, for completeness a compact form of the equations used for the teetering rotor are reproduced here:
A1and
A 2where
where *N*_{ψ,i} and *M*_{β,i} are the overall torque and flapping moment, respectively, on blade *i* and *N*_{aero,i} and *M*_{aero,i} are the aerodynamic torque and flapping moment, respectively, for each blade *i*. *M*_{blade} and *I*_{blade} are the blade mass moment and second moment of inertia, respectively. , and are the angular velocities of the rotor hub in the shaft axes. *a*^{shaft}_{hub,i} and *α*^{shaft}_{hub,i} are the translational and angular acceleration vectors of the hub in the shaft axes, respectively. *α*^{disc}_{hub} is the angular acceleration vector of the hub in the rotor disc axes. Finally, the state values are allocated for each rotor blade according to table 1. Note that moments at the hub about either the teetering or rotation axes due to factors such as friction are assumed negligible at this stage.

A two-dimensional individual blade element approach is adopted. Each blade is divided into an *N* number of elements and then the aerodynamic forces for each element are calculated numerically. The lift and drag forces acting on each blade element can be computed by first evaluating the elemental flow velocities and angles of attack (figure 12). The tangential and perpendicular components of resultant velocity at a blade element are given by
A 3and
A 4where *V*_{elem} and *W*_{elem} are the translational velocities of the blade element in the rotation and flap directions, respectively, owing to the kinematic motion of the blade. *ν*_{i} is the induced velocity at the blade element. The local angle of attack is given by
A 5where *θ*_{elem} is the elemental pitch angle and *ϕ*_{elem} is the elemental inflow angle, which can be calculated as follows:
A 6where atan2 is the ‘quadrant-arctangent’ function (atan2 function in Matlab). The local pitch angle *θ*_{elem} is a combination of all the local pitch angle contributions, in this case
A 7where *θ*_{col} and *θ*_{pre.twist,elem} are the blade collective pitch and the blade local built-in twist, respectively. The elemental lift, drag forces and pitching moment are then calculated assuming two-dimensional steady flow.
A 8and
A 9where *ρ*, *c*_{elem}, , *δr*_{elem}, *C*_{L,elem} and *C*_{D,elem} are the local air density, the blade elemental chord, the elemental resultant flow velocity, the width of the blade element, and the lift and drag coefficients at each element, respectively. *α*_{elem} and *Re*_{elem} are the elemental angle of attack and Reynolds number, respectively. The aerodynamic loads and moments for each individual blade element are calculated numerically using nonlinear look-up tables for lift and drag coefficients. Experimental data for a NACA0015 aerofoil section are used for each blade element [42]. The local forces acting on a blade element in the blade coordinate system can therefore be determined from the elemental lift and drag forces.
A 10

Finally, the aerodynamic forces and moments acting upon the whole blade at the hub are calculated by summing up the whole elemental forces and moments for each blade.
A 11and
A 12where is the blade element position vector relative to the rotor hub. The thrust *T*_{aero} and lift *L*_{rotor} generated by the rotor can be calculated as follows:
A 13and
A 14

The three-state inflow model is given as
A 15where *ν*_{i} is the induced velocity at the *i*th element of radius *r* and *R* is the blade tip radius. The induced velocity components *ν*_{0}, *ν*_{s} and *ν*_{c} are given in the wind axes by
A 16where
A 17and
A 18*T*_{aero}, *L*_{aero} and *M*_{aero} are the thrust and the aerodynamic lateral and pitching moments, respectively, in the wind axes and the skew angle *χ* can be obtained from
A 19where the ‘quadrant-arctangent’ function (atan2 function in Matlab) is used. *ν*_{T} and *ν*_{m} are expressed as
A 20and
A 21where *u*_{hub}, *v*_{hub} and *w*_{hub} are the flow velocity components relative to the rotor hub, and *ν*_{i,m} is the induced velocity from momentum theory, which can be evaluated using the Glauert equation [51].

## Footnotes

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

- Accepted July 8, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.