2.1 9-DoF whirl flutter model

An 18-state, 9-DoF model originally presented by Johnson [Reference Johnson25] is used. A schematic of the model to be used is shown in Fig. 2. A rotor of radius R spins at angular velocity Ω, attached to a shaft via a gimballed hub that allows it to rotate in pitch and yaw, forming respective angles $\beta{\kern-1pt}$ _1C and $\beta{\kern-1pt}$ _1S between the hub plane’s instantaneous position and undeflected position. The shaft of length h, is also able to move within the reference frame coordinate system (x, y, z) as indicated, by virtue of movement of the supporting wing structure in flapwise bending q ₁ , chordwise bending q ₂ and torsion p. The wing is cantilevered about its root, and has span y _Tw and chord c _w and in its undeflected state the rotor shaft points in the +z direction. A freestream velocity V inward along – z is incident on the rotor, and aggregate stiffness and damping properties exist at both the gimbal and the wing. Additionally, the rotor blades are able to move in cyclic lead-lag, causing rectilinear motion of the CG in both the lateral ( $\zeta$ _1C ) and vertical ( $\zeta$ _1S ) directions as viewed in the x-y plane. The rotor blades may also flap collectively, constituting a coning motion $\beta{\kern-1pt}$ ₀ , and lead-lag collectively, constituting a perturbation of rotor speed $\zeta$ ₀ . The aerodynamics of both the blades and the wing are modelled using quasi-steady strip theory [Reference Ribner26]. Important to note is that all blade deflection is assumed to be rigid body rotation about the hub. That is, this motion is a simplified aggregation of deflections that might arise from both hub gimballing and from elastic deformation of the blades’ flex-beams. Typically these motions are understood to be separate in real-world tiltrotor systems and models made of them. Only the first mode of the blade motions (both flap and lead-lag) and the aforementioned wing motions are considered, due to the assumption that higher modes have negligible participation in the coupled wing-rotor motion. Additionally, the motions are considered to be uncoupled to each other. The equations of motion of the model are far too long to state in full here, though they are of the form:

(1) \begin{equation}{{\bf{M}}\ddot{\textbf{X}}} + {\bf{C}}\ddot{\textbf{X}} + {\bf{KX}} = \textbf{0}\end{equation}

where M is the mass matrix, C is the damping matrix and K is the stiffness matrix. The latter matrices contain both structural and aerodynamic terms. X is the vector of generalised displacements, containing the various system-defining quantities discussed earlier:

(2) \begin{equation}{\bf{X}} = {\left[ {{\beta _{1C}}\quad{\beta _{1C}}\quad{\zeta _{1C}}\quad{\zeta _{1C}}\quad{\beta _0}\quad{\zeta _0}\quad{q_1}\quad{q_2}\quad p} \right]^T}\end{equation}

The parameter values used in the model are shown in Table 1. The nondimensionalisation factor is generally $I_{b}\Omega^{2}$ ; details of this and the few exceptions are given in Ref. Reference Johnson25.

Figure 2. Schematic diagram of gimballed rotor-nacelle model.

Table 1. Datum parameter values used in the model

2.2 Nonlinear adaptation

In the original model’s derivation [Reference Johnson25], within the rotating frame, the rotor blades’ (linear) flapping and lead-lag stiffnesses are specified implicitly in terms of the per-rev frequencies of each uncoupled motion, $\nu_{\beta}$ and $\nu_{\zeta}$ respectively, which the model specifies as functions of μ and Ω to account for centrifugal stiffening and aerodynamic effects. The elastic structural contributions are therefore $\nu^{{2}}_{\beta}{\beta_{{\textit{m}}}}$ in the flapping equations and $\nu^{{2}}_{\zeta}{\zeta_{{\textit{m}}}}$ in the lead-lag equations, where ${\beta_{{\textit{m}}}}$ is the instantaneous flap angle of the m ^th blade and ${\zeta _m}$ is the instantaneous lead-lag angle of the m ^th blade. As mentioned before, the flapping and lead-lag motions are modelled as rigid-body aggregations of motion, defined by these elastic properties.

In this research, the aforementioned linear elastic structural contribution in the flapping equations is replaced with a cubic polynomial of the form:

(3) \begin{equation}\;{K_1}{\beta _m} + {K_2}\beta _m^3\end{equation}

where the first term is simply the linear expression used in the original model. That is, K ₁ is equal to $\nu^{{2}}_{\beta}$ . This nonlinear stiffness expression can provide hardening behaviour by using a positive value of K ₂ , and softening behaviour by using a negative value. In this paper, a positive value of 10 was used to allow a focus on hardening behaviours, though the use of a softening profile with this particular rotor-nacelle model was explored extensively in Ref. Reference Mair, Rezgui and Titurus24.

The nonlinear stiffness profile compared to the original linear model (i.e. K ₂ = 0), is demonstrated in Fig. 3.

Figure 3. Nonlinear stiffness profile (green) compared to linear profile (blue).

To obtain the dimensional values of K _1,2 they should be multiplied by I _b $\Omega$ ² (see Table 1). K ₁ is selected as the independent variable in the results to follow and it is therefore the continuation parameter in the continuations shown. However, it will be controlled not by its absolute value but rather as a multiple of its datum value, $\nu _{\beta ,datum}^2$ . This variable, termed $K^{\#}_{\beta}$ , is defined as $\frac{{{K_1}}}{{\nu _{\beta ,datum}^2}}$ . Lastly, as this quantity is controlled indirectly by adjusting the frequency $\nu_\beta$ in the model, the variable $K^{\#}_{\beta}$ is referred to as the effective blade flapping stiffness.

In Ref. Reference Johnson25 the rotor equations of motion are obtained by considering flapping and lead-lag equations for each individual blade within the rotating frame of reference, and summing these over the whole rotor via a Fourier Coordinate Transform to obtain equations in the external, non-rotating frame. These summations also lead to a transformation of the degrees of freedom: coning angle $\beta{\kern-1pt}$ ₀ , gimbal pitch $\beta{\kern-1pt}$ _1C and gimbal yaw $\beta{\kern-1pt}$ _1S . These are defined thus:

(4) \begin{equation}\;{\beta _0} = \frac{1}{N}\mathop \sum \limits_{m = 1}^N {\beta _m}\end{equation}

(5) \begin{equation}\;{\beta _{1C}} = \frac{2}{N}\mathop \sum \limits_{m = 1}^N {\beta _m}\cos {\psi _m}\end{equation}

(6) \begin{equation}\;{\beta _{1S}} = \frac{2}{N}\mathop \sum \limits_{m = 1}^N {\beta _m}\sin {\psi _m}\end{equation}

where $\psi$ _m is the instantaneous angle of the m ^th blade. That is, by applying the respective summation in each of Equations (4)–(6) to the blade flapping equations in the rotating frame, the equations of motion in the non-rotating frame are obtained. A similar process is used to form the lead-lag equations of motion in the non-rotating frame.

As the $K_{1}\beta_{m}$ term of the new nonlinear profile is just the existing $\nu _\beta ^2$ term in the flapping equations, only the new K ₂ $\beta _m^3$ term needs to be processed in this way. In a more sophisticated model, the blade flapping and coning motions would have differing origins, but due to the aggregated rigid body modelling used in Johnson’s model, the stiffness nonlinearity introduced here will manifest in both the gimballing and coning equations. Using the following relationship:

(7) \begin{equation}\;{\beta _m} = {\beta _0} + {\beta _{1C}}\cos {\psi _m} + {\beta _{1S}}\sin {\psi _m}\end{equation}

the K ₂ $\beta _m^3$ term can be expanded:

(8) \begin{align}{K_2}\beta _m^3 & = {K_2}{\left( {{\beta _0} + {\beta _{1C}}\cos {\psi _m} + {\beta _{1S}}\sin {\psi _m}} \right)^3}\nonumber\\[4pt]& = {K_2}\left[\beta _0^3 + 3\beta _0^2{\beta _{1C}}\cos {\psi _m} + 3\beta _0^2{\beta _{1S}}\sin {\psi _m} + 3{\beta _0}\beta _{1C}^2{\cos ^2}{\psi _m} + 6{\beta _0}{\beta _{1C}}\cos {\psi _m}{\beta _{1S}}\sin {\psi _m}\right.\\[4pt] & \quad + 3{\beta _0}\beta _{1S}^2{\sin ^2}{\psi _m} + \beta _{1C}^3{\cos ^3}{\psi _m} + 3\beta _{1C}^2{\cos ^2}{\psi _m}{\beta _{1S}}\sin {\psi _m} + 3{\beta _{1C}}\cos {\psi _m}\beta _{1S}^2{\sin ^2}{\psi _m}\nonumber\\[4pt]& \quad \left. + \beta _{1S}^3{\sin ^3}{\psi _m}\right]\end{align}

The application of the respective summation in each of Equations (4)–(6) produces the additions to be made to the coning ( ${\beta _0}$ ) equation, the gimbal pitch ( ${\beta _{1C}}$ ) equation and the gimbal yaw ( ${\beta _{1S}}$ ) equation. The respective additions to these equations are:

(9) \begin{equation}\frac{{{K_2}}}{4}\left( {2\beta _0^2 + 3\beta _{1C}^2 + 3\beta _{1S}^2} \right){\beta _0}\end{equation}

(10) \begin{equation}\frac{{{K_2}}}{4}\left( {12\beta _0^2 + 3\beta _{1C}^2 + 3\beta _{1S}^2} \right){\beta _{1C}}\end{equation}

(11) \begin{equation}\frac{{{K_2}}}{4}\left( {12\beta _0^2 + 3\beta _{1C}^2 + 3\beta _{1S}^2} \right){\beta _{1S}}\end{equation}

The model in both its linear and nonlinear forms was implemented in MATLAB R2015a [27].

3.1 Eigenvalue analysis

To form a baseline stability analysis of the model before the nonlinear analysis is conducted, eigenvalue analysis is employed. This standard method for analysing linear dynamical systems uses the Jacobian matrix ${\bf{J}}$ , defined as:

(12) \begin{equation}{\dot{\bf{Y}}} = {\bf{JY}}\end{equation}

where Y, the state vector, is defined as:

(13) \begin{equation}{\dot{\bf{Y}}} = \boldsymbol{{f}}\left( {{\bf{Y}},{\bf{p}}} \right),{\rm{\;}}{\bf{Y}} \in {\mathbb{R}^{2n}},{\rm{\;}}{\bf{p}} \in {\mathbb{R}^l}\end{equation}

(14) \begin{equation}{\bf{Y}} = \left[ {\begin{array}{*{20}{c}}{\bf{X}}\\[3pt]{\dot{\bf{X}}}\end{array}} \right]\end{equation}

where $\boldsymbol{{p}}$ is a vector of $\boldsymbol{{l}}$ system parameters and $\boldsymbol{{n}}$ is the dimension of ${\bf{X}}$ , the vector of generalised displacements defined in Equation (2).

The $2n$ eigenvalues of the Jacobian matrix contain information about the decay rate (i.e. stability) and frequency of the system’s modes. The undamped natural frequency $\omega $ and damping ratio $\zeta $ for each mode are calculated from the real and imaginary parts of its eigenvalue $\lambda $ using Equations (15) and (16).

(15) \begin{equation}\omega = \sqrt {{\mathop{\rm Re}\nolimits} {{\left( \lambda \right)}^2} + {\rm{Im}}{{\left( \lambda \right)}^2}} \end{equation}

(16) \begin{equation}\zeta = \frac{{ - {\mathop{\rm Re}\nolimits} \left( \lambda \right)}}{\omega }\end{equation}

Scripts for this eigenvalue analysis were written in MATLAB so that a direct interface with the model was possible.

3.2 Continuation and bifurcation methods

For nonlinear systems, numerical continuation and bifurcation theory are used. Continuation is a numerical method that finds the steady-state solution values of a dynamical system as one or more of its parameters, called the continuation parameter(s), is varied [Reference Rezgui, Lowenberg, Jones and Monteggia28], constructing solution branches or “continuing” the set of solutions. Continuation can find both fixed point solutions and periodic solutions. At fixed points the system is considered to be in equilibrium, such as a rigid pendulum standing motionless at either the bottom or top of its arc of motion. Periodic solutions, also known as limit cycle oscillations (LCOs), are closed trajectories through the state space that return precisely to their starting point and constitute motions that repeat periodically. For each solution point calculated, the stability is then computed. For fixed points, an eigenvalue analysis of the type described in Section 3.1 can be used, requiring local linearisation in the case of a nonlinear system. Periodic behaviour on the other hand requires Floquet theory to determine stability [Reference Kuchment29].

A bifurcation is a qualitative change in the system behaviour due to the variation of a parameter. When the stability of a system changes, or the type of the solution changes (fixed/periodic), the system bifurcates. The points at which these stability changes happen are called bifurcations. If the system is nonlinear, new solution branches may emerge from the bifurcation points, leading to the presence of multiple solutions for a given set of system parameters.

Branches are explored as desired to uncover the global dynamics of the system within the desired parametric ranges. The results of the continuations are displayed on bifurcation diagrams, where the values of solution branches are shown as the continuation parameter value varies. The type (equilibrium/periodic) of each solution branch, along with the location of any bifurcations it encounters, is also indicated. The solutions exist in a space whose number of dimensions is the number of states plus the number of continuation parameters. To allow plotting, conventional bifurcation diagrams are 2D graphs with the continuation parameter on the x-axis and a chosen state on the y-axis. The plotting of the solutions in terms of that chosen state is known as a “projection” or a “plane”, e.g. the “ ${\beta _{{1}{\textit{C}}}}$ projection”, or the “ ${{\textit{K}}_{{\textit{q}2}}}\hbox{-}{\textit{p}}$ plane”. Alternatively, if a 2-parameter continuation is conducted, the results can be plotted on a bifurcation diagram where both axes are parameters.

Bifurcation diagrams for this research were produced using the Dynamical Systems Toolbox for MATLAB by Coetzee [Reference Coetzee, Krauskopf and Lowenberg30], which uses an implementation of AUTO-07P [Reference Doedel, Fairgrieve and Sandstede31]. Time simulations were also used to corroborate the predictions of both stability methods.

Key bifurcation types that are relevant to understanding the behaviour of this rotor-nacelle system, are Hopf bifurcations and torus bifurcations. At a Hopf bifurcation, the stability of a fixed point (i.e. an equilibrium) changes, and a periodic solution arises, caused by a pair of complex conjugate eigenvalues crossing the complex plane imaginary axis. At a torus bifurcation, a torus-shaped manifold smoothly emerges around an existing periodic solution, encasing it. While the periodic solution still exists, trajectories may now flow on the torus manifold and crucially, such trajectories are quasi-periodic. That is, the motion around the torus never quite exactly repeats. Both these bifurcation types are illustrated and explained in the context of this research in Section 4. A generalised, in-depth theoretical treatment of these bifurcations is given in Ref. Reference Thompson and Stewart32.

4.1 Linear analysis

The eigenvalue analysis methods detailed in Section 3.1 are now applied to the model to fulfil two purposes. Firstly, this work’s implementation of Johnson’s model is validated by conducting a parameter sweep identical to one shown in Johnson’s original work [Reference Johnson25]. Secondly, the eigenvalue analysis is used to generate a baseline stability boundary that will be the basis of comparison for understanding the effects of the presence of the nonlinearity introduced.

Shown in Fig. 4 are the eigenvalues and associated modal damping ratios and modal frequencies for the system undergoing an airspeed sweep of V $\in$ [13, 309]m/s. All other parameters are kept at their datum values, as given in Table 1. This figure corresponds to Fig. 18 within Ref. Reference Johnson25, which uses an airspeed range of 25–600Kn (i.e. 13–309m/s), uses the damped definition of natural frequency (i.e. $\omega = {\rm{Im}}(\lambda)$ rather than Equation (15)), only shows the positive-imaginary root of each complex conjugate pair, and uses two co-located x-axes of different scales. The plot lines in this original figure have been imported to Fig. 4 and are shown in thin, beaded black lines. In the original figure, the damping ratio data is only provided in the bracket of ζ $\in$ [0,0.15]. Furthermore, the eigenvalue loci data for the complete sweep are not shown in all modes. All the data that is given is shown here.

Figure 4. Eigenvalues (top) and their corresponding modal damping ratios (centre) and frequencies (bottom) for a sweep in airspeed V, original linear model.

The mode naming used in Fig. 4 (shown far right) is taken directly from Johnson, who devised the naming scheme based on the prominence of participation of the system’s degrees of freedom in each mode, and the proximity of modal frequencies to uncoupled natural frequencies of the system’s degrees of freedom.

The concept of a stability boundary diagram between two parameters can be useful for understanding the sensitivity of a system’s stability to changes in various parameters. Such a diagram can be produced from a grid of the combinations of different values for each parameter. The Jacobian matrix and its eigenvalues are calculated at each point, and a surface is overlaid on the grid whose height at each point is determined by the maximum real component of the Jacobian’s eigenvalues there. As only one unstable eigenvalue is required for overall system instability, a horizontal plane cut of this surface at the level 0 will produce a contour that denotes the boundary between the stable and unstable regions of the grid.

As aeroelasticity is such an important consideration during the design of tiltrotor aircraft, stability boundaries between structural parameters that are controllable to some degree are a useful tool for the designer. To provide a simple basis for the introduction of the continuation methods employed in subsequent sections, the stability boundary chosen for this work is between (linear) gimbal flapping stiffness $K^{\#}_{\beta}$ (which is used as the continuation parameter in Section 4.2) and the datum-normalised wing torsional stiffness $K_p^\# = \frac{{{K_p}}}{{{K_{p,datum}}}}$ . The latter is chosen as the boundary thereby lends insight into the relative influence of the two main components of the model: the rotor and the wing. This stability boundary, generated with the linear analysis detailed in Section 3.1, is shown in Fig. 5. All other system parameters use their datum values, as given in Table 1.

Figure 5. Stability boundary between datum-normalised blade flapping stiffness $K^{\#}_{\beta}$ and datum-normalised wing torsional stiffness $K^{\#}_{p}$ .

Crucially, the linear analysis returns the same results for both the original linear form of the model and the adapted nonlinear version, regardless of the value of K ₂ used. For the fixed point solution branches, this is correct. However, as will be shown in Section 4.3, nonlinearity is able to create new periodic solutions whose parametric ranges may not match those of the fixed point branches. The linear analysis is not able to detect the new periodic solutions, and therefore it is not sufficient just to use a model that successfully replicates nonlinearities in a real-world system: with only linear analysis, the impacts of these nonlinearities are still lost.

4.2 Bifurcation analysis: Linear model

To demonstrate the link between eigenvalue analysis and CBM, Fig. 5 can also be generated by continuation methods, as the system has an equilibrium at Y = 0 that can be used as a starting solution. This vector corresponds to the rotor-nacelle system at rest and in an entirely undeflected state. Generating the stability boundary using CBM in fact affords deeper insight than the contour cut method.

Firstly, a one-parameter continuation is performed on the system. Choosing $K^{\#}_{p}$ = 0.55 so that a continuation in $K^{\#}_{\beta}$ will intersect a region of interest in the stability boundary, the bifurcation diagram shown in Fig. 6 (top) is obtained. The modal damping ratios for this sweep are shown at the bottom. A key to the bifurcation diagrams in this paper is shown in Table 2.

Figure 6. Bifurcation diagram (top), $\beta$ _1C projection, with corresponding modal damping ratios (bottom), as $K^{\#}_{\beta}$ is varied while $K^{\#}_{p}$ = 0.55.

Table 2. Key to symbols and line colours used in bifurcation diagrams in this paper

Moving in the sense of decreasing blade flap stiffness on the damping ratio plot, the q ₁ mode (red) is the first to lose stability, followed by the p mode as previously mentioned, followed by the q ₂ mode (yellow). However the q ₁ mode restabilises shortly afterward, and by approximately $K^{\#}_{\beta}$ = 0.45 the q ₂ mode has also restabilised. While this variety of types of whirl flutter is not in itself a hazard, it is their coexistence in the parameter space that may pose a threat, as if more than one stable whirl flutter solution exists at a certain value, it increases the likelihood that a perturbation will cause the system to encounter whirl flutter of some form. The continuation finds five Hopf bifurcations (hollow square icons in the top of the figure) and when their $K^{\#}_{\beta}$ coordinates are cross-referenced with the damping ratio plot, it becomes apparent that they may signal either the loss of stability or the re-stabilisation of an individual whirl flutter mode. The continuation also finds that the gimbal pitch coordinate ( $\beta{\kern-1pt}$ _1C ) of this equilibrium remains at 0° over the whole continuation, however only the branch segment to the right of the right-most Hopf bifurcation at approximately $K^{\#}_{\beta}$ = 0.65 is stable. This correlates with the damping ratio plot, which shows that there is always at least one unstable eigenvalue at $K^{\#}_{\beta}$ < 0.65. The equilibrium remaining at 0° is to be expected given the linear form of the equations.

Performing two-parameter continuations in $K^{\#}_{\beta}$ and $K^{\#}_{p}$ on each of these Hopf bifurcations reconstructs the stability boundary and this is shown in Fig. 7. The very bottom-right of the boundary is defined by a branch point bifurcation of the pitchfork type that is found by a separate continuation downward in $K^{\#}_{p}$ at any fixed value of $K^{\#}_{\beta}$ (not shown). This pitchfork gives rise to secondary equilibrium branches though as they do not relate to whirl flutter, they are not explored. Regions shaded red are subject to oscillatory instability (an unstable complex conjugate root pair) and regions shaded blue are subject to non-oscillatory instability (an unstable single real root).

Figure 7. $K^{\#}_{\beta}$ - $K^{\#}_{p}$ stability boundary regenerated through two-parameter continuations.

Two levels of $K^{\#}_{p}$ are chosen as analysis cases for the rest of the chapter: $K^{\#}_{p}$ = 1.1 and $K^{\#}_{p}$ = 0.2. They are also indicated in Fig. 7. This choice of cases provides a concise demonstration of the influence of $K^{\#}_{p}$ as they span the range of analysis: one at high pitch stiffness and one at lower stiffness.

4.3 Bifurcation analysis: Nonlinear model

Using the nonlinear stiffness profile given in Equation (3), with a K ₂ value of 10, a continuation in $K^{\#}_{\beta}$ for the first analysis case ( $K^{\#}_{p}$ = 1.1) is conducted and shown in Fig. 8. Two projections are shown to illustrate the behaviour in the two main components of the model, the rotor (gimbal pitch $\beta{\kern-1pt}$ _1C , top) and the wing (wing torsion p, bottom).

Figure 8. Bifurcation diagram for case 1 ( $K^{\#}_{p}$ = 1.1), $\beta$ _1C and p projections.

It is worth noting that several of the solutions found here are likely to have too large amplitude to be physically meaningful. That is, a real-world system would experience structural failure before achieving such displacements. As the model does not consider structural strength, such failure may only be implied or assumed. However, they are still useful for clear demonstration of the key phenomena described in this section, which are likely to occur in smaller-amplitude solutions in either this model or similar systems.

Once again, an equilibrium branch that remains at 0° is present that is visible in both projections. This branch will hereafter be referred to as “the main branch”. Three Hopf bifurcations are visible on this main branch, whose locations in $K^{\#}_{\beta}$ correspond to the three intersections of a line at $K^{\#}_{p}$ = 1.1 with the Hopf loci present at that level (see Fig. 7).

New to this bifurcation diagram, compared to that shown in Fig. 6, is the indication of the periodic solution branches that emanate from the Hopf bifurcations. These are indicated slightly differently to equilibrium branches; while an equilibrium has only one coordinate in each state at each continuation parameter value, the same is not true for periodic solution branches as at each parameter point, a whole LCO exists, covering a range of values in the states. The convention then is to indicate periodic solution branches via the maximum (i.e. most positive) value of the given projection state in the LCO at each point.

The most significant wing torsion participation occurs in the periodic solution branch emanating from the middle Hopf bifurcation, with much less prominent participation in the branches emanating from the left and right Hopf bifurcations. The hardening component in the blade flap stiffness profile bends the periodic solutions leftward. As $K^{\#}_{\beta}$ is lowered, the amplitude of all three whirl flutter branches increases in both projections. Lastly, there are three torus bifurcations present on the diagrams. As explained earlier, torus bifurcations signal the emergence of a torus manifold from an LCO, on which system trajectories may flow in a quasi-periodic manner. This is explored and illustrated shortly.

The periodic solution branches mostly coexist – that is, at a given continuation parameter value, more than one periodic solution branch is generally present – and the system may join any stable solution branch segments depending on what perturbation it receives. This constitutes the possibility of whirl flutter. It is important to note that CBM’s prediction of stable branches in parametric regions that lie within the unstable region predicted by linear analysis does not indicate that these regions have been made safe by the presence of the nonlinearity. The model does not account for any damage that the occurrence of whirl flutter LCOs would likely cause to the system, leading to the widening of the oscillations and likely structural failure.

Bifurcation diagrams, while not counterintuitive, require some familiarisation to be used effectively. Thus, to aid the reader in becoming fluent at visualising Fig. 8 and the remaining bifurcation diagrams in this paper, a selection of $\beta{\kern-1pt}$ _1S - $\beta{\kern-1pt}$ _1C phase planes is provided in Fig. 9, showing a selection of the various types of solutions, both stable and unstable, present at a number of $K^{\#}_{\beta}$ cuts.

Figure 9. Phase plane planes showing solutions in detail at various $K^{\#}_{\beta}$ values in Fig. 8.

Observing the figure from right to left, the first cut at $K^{\#}_{\beta}$ = 0.8 shows the stable main branch with no other solutions present, the second at $K^{\#}_{\beta}$ = 0.45 shows the main branch now unstable with a single stable whirl flutter LCO present, and the third shows three whirl flutter LCOs all present at the $K^{\#}_{\beta}$ value of 0.1. Some torus flow was located at this value of $K^{\#}_{\beta}$ on the branch emanating from the right-most Hopf bifurcation at $K^{\#}_{\beta}$ = 0.497, shortly after the torus bifurcation has made the LCO unstable. The LCO at this point is visible on the left-most phase plane of Fig. 9; it is the larger of the two unstable LCOs (red). This torus flow is demonstrated in Fig. 10 in a 3-dimensional projection of the phase space using gimbal pitch $\beta{\kern-1pt}$ _1C , gimbal yaw $\beta{\kern-1pt}$ _1S and wing torsion p, produced by a time simulation over a period of 250 rotor revolutions (bottom). A Poincaré section placed at $\beta{\kern-1pt}$ _1C = 0° on the positive $\beta{\kern-1pt}$ _1S side of the LCO shows a cross-section of the torus’ 3D projection at this point in the LCO, over a period of 5,000 rotor revolutions (top). The torus can clearly be seen surrounding the LCO.

Figure 10. Torus flow found at $K^{\#}_{\beta}$ = 0.1 in case 1, shown in ( $\beta$ _1C, $\beta$ _1S, p) space (bottom), with Poincaré section placed at $\beta$ _1C = 0° (top).

The second analysis case ( $K^{\#}_{p}$ = 0.2) is now considered, shown in Fig. 11. Once again, there are three Hopf bifurcations present. The middle and left periodic solution branches are both unstable within the domain of analysis. The right-most Hopf’s flutter branch is the only one of the three to contain stable portions, and after emanating from the main branch it folds back and forth in $K^{\#}_{\beta}$ a number of times before eventually departing to the left. The much lower wing torsional stiffness value of case 2 allows all three whirl flutter branches to be much larger in amplitude than those in case 1.

Figure 11. Bifurcation diagram for case 2 ( $K^{\#}_{p}$ = 0.2), $\beta$ _1C and p projections. A time simulation starting close to the LCO branch at $K^{\#}_{\beta}$ = 1.015 is also shown.

There are also two torus bifurcations present, and the time history of the $\beta{\kern-1pt}$ _1C state during some torus flow that was found is shown at the top of Fig. 11. The quasi-periodicity of the second half of the time history is evident from the notable breaking of the strict regularity symptomatic of LCO activity that is present in the first half of the time history.

Most pressing however is the overhanging of the stable main (equilibrium) branch by a portion of the first Hopf’s whirl flutter branch which, due to the torus bifurcation at approximately $K^{\#}_{\beta}$ = 1, is surrounded by torus flow in at least some neighbourhood of that torus bifurcation. A time simulation with $K^{\#}_{\beta}$ = 1.015 (top of Fig. 11) is shown to demonstrate the system’s ability to be attracted to whirl flutter behaviours – of both the LCO and quasi-periodic varieties – in this vicinity. The simulation of 1,000 rotor revolutions is started with initial conditions 99% the size of a point on the LCO (in all states) and the transition from the LCO to the quasi-periodic behaviour on the torus surrounding the LCO is noticeable in the time histories for each state shown. This overhang stretches as far in $K^{\#}_{\beta}$ as the torus bifurcation, and is partially shadowed by a lesser overhang by the same branch at a larger amplitude, which provides stable LCOs up to approximately $K^{\#}_{\beta}$ = 0.95.

4.4 Impacts on stability boundary

The whirl flutter branches created by the nonlinearities are seen to overhang the stable main branch (the undeflected position of the rotor-nacelle system) at several points in the results. This means that flutter can be encountered in parametric regions that linear eigenvalue analysis predicts to be stable. The stability boundary must therefore be redrawn to include such regions.

To do so, individual instances of overhang observed in the bifurcation diagrams produced for the two $K^{\#}_{p}$ analysis cases must be tracked over the full $K^{\#}_{\beta}$ - $K^{\#}_{p}$ domain under consideration. This is done through a variety of different continuations: continuations in $K^{\#}_{\beta}$ at levels of $K^{\#}_{p}$ that are otherwise not used for analysis cases, continuations in $K^{\#}_{p}$ at various positions in $K^{\#}_{\beta}$ , and two-parameter continuations in $K^{\#}_{\beta}$ and $K^{\#}_{p}$ simultaneously.

There are two instances of overhang that require tracking, shown in Fig. 12, top. The first is an LCO overhang, marked ‘O1’. In case 1 it originates from the middle Hopf bifurcation (the connection is not visible within the $K^{\#}_{\beta}$ range shown), and overhangs as far as approximately $K^{\#}_{p}$ = 0.76, while in case 2 it originates from the right-most Hopf bifurcation and overhangs as far as approximately $K^{\#}_{p}$ = 1. The other overhang present is an outcrop of the rightmost whirl flutter branch, which contains a torus bifurcation, marked ‘O2’ in Fig. 12. It does not have a corresponding presence in the previously shown case 1 diagrams due to not having any Hopf connections there, and exists there solely due to a large overhang in $K^{\#}_{p}$ . Tracking these overhangs within the $K^{\#}_{\beta}$ - $K^{\#}_{p}$ plane renders the redrawn stability diagram shown at the bottom of Fig. 12.

Figure 12. Overhang instances (top) with redrawn $K^{\#}_{\beta}$ - $K^{\#}_{p}$ stability boundary (bottom).

The new hatched regions drawn on Fig. 12 represent parametric regions in which whirl flutter is now possible due to the presence of nonlinearities. Though they occupy what was termed the “unstable” region on the original stability boundary shown in Fig. 5, these new hatched regions cannot strictly be termed “unstable” themselves as the whirl flutter hazard they constitute is caused by the presence of stable LCOs. Instead, the more accurate term “unsafe” is used.

4.5 Analysis of dynamical behaviours

The incursions into the supposedly stable region of the stability boundary are qualitatively different: one is an LCO and the other is a torus flow. While the amplitude of the solutions found in this paper is too large for this qualitative difference to be meaningful, smaller solutions found in similar systems could benefit from insight that can be gained here into the implications of this difference. Some further analysis can be employed to understand the differences between these two dynamical behaviours with a view to understanding their respective whirl flutter implications. The first method, the use of a Poincaré section, was introduced in Section 4.3. The second is the use of the Fourier transform, using MATLAB’s inbuilt fft function, to provide information on the frequency components of the two behaviours.

4.5.1 LCO analysis

An LCO point is analysed first, shown in Fig. 13. The chosen point is on the O1 overhang in case 2, at approximately $K^{\#}_{\beta}$ = 0.95. The LCO (dashed red), as calculated by AUTO, along with a time simulation left to run in its near vicinity (blue), are shown in ( $\beta{\kern-1pt}$ _1C , $\beta{\kern-1pt}$ _1S , p) space. The simulation follows the LCO perfectly, passing around it a number of times. This is reflected in the Poincaré section (bottom left), where only a single repeated blue dot is visible, at the exact intersection point of the LCO path. The frequency spectrum shows a number of regularly spaced peaks that are consistent with the use of a nonlinear stiffness profile [Reference Thompson and Stewart32]. These harmonics are directly related to the use of a cubic polynomial.

Figure 13. Analysis of an LCO found in the system. Phase space (top), Poincaré section (bottom left) and frequency analysis (bottom right).

4.5.2 Quasi-periodic analysis

By comparison, the time simulation left to run in the vicinity of the LCO shown in Fig. 14 finds the torus manifold that has emanated from the torus bifurcation shown in the O2 overhang, also in case 2. The frequency spectrum is essentially broadband white noise with a few peaks superimposed upon it.

Figure 14. Analysis of a torus flow found in the system. Phase space (top), Poincaré section (bottom left) and frequency analysis (bottom right).

The quasi-periodic motion is significantly more complex than the LCO and has a wide range of constituent frequencies. Compared to the LCO from which a torus flow originates, the torus flow also has larger amplitude. These two aspects potentially make the quasi-periodic motion more of a whirl flutter hazard than an LCO.