Stability and dynamical analysis of whirl flutter in a gimballed rotor-nacelle system with a smooth nonlinearity

Abstract

Whirl flutter is an aeroelastic instability that affects aircraft with propellers/rotors. With their long and flexible rotor blades, tiltrotor aircraft are particularly susceptible. Whirl flutter is known to have destroyed aircraft and in the best case it constitutes a fatigue hazard. The complexity of whirl flutter analysis increases significantly with the addition of nonlinearities, due to the more complex dynamical behaviours that emerge as a result. Most whirl flutter stability analyses in current literature are grounded in linear theory, preventing the full discovery of the nonlinearities’ effects. Continuation and bifurcation methods (CBM) may instead be used to fully appreciate and analyse the effects of the presence of nonlinearities. Previous CBM-based work on nonlinear gimballed hub rotor-nacelle models, representing those found on tiltrotor aircraft, are capable of whirl flutter in parametric regions declared safe by linear analysis. Furthermore, it was found that they are capable of complex behaviours including limit cycle oscillations, quasi-periodic behaviour and even chaos, though the whirl flutter implications of such behaviours has not been explored. This paper investigates the impact of a smooth structural nonlinearity on the whirl flutter stability of a basic gimballed rotor-nacelle model, compared to its baseline linear stiffness version. A 9-DoF model with quasi-steady aerodynamics, a flexible wing and blades that can move both cyclically and collectively in both flapping and lead-lag motions, producing gimbal flap-like behaviour, was adopted from existing literature. A smooth stiffness nonlinearity was introduced in the blade flapping stiffness and CBM was used to find the new whirl flutter behaviours created by the presence of the nonlinearity. Time simulations, Poincaré sections and spectral analysis were then used to investigate the various behaviours found. This in turn allowed recommendations to be made concerning preferable and/or hazardous parameter combinations of use to the tiltrotor designer.

Abbreviations

DoF

degree(s) of freedom

LCO

limit cycle oscillation

R, [rad/s]

Rotor angular velocity

V, [m/s]

Freestream velocity

$\gamma$ , [-]

Lock number

$\mu$ , [-]

Advance ratio

K, [-]

Normalised stiffness

C, [-]

Normalised damping

$I^{(\ast)}$ , [-]

Normalised inertia

$\nu$ , [-]

Per-rev frequency

1.0 Introduction

Tiltrotor aircraft, such as the XV-15, shown in Fig. 1, are growing in importance due to their unique flight envelope. A principal tiltrotor design challenge is ensuring that the airframe structure prevents flutter from occurring within the operating envelope. The most basic form of whirl flutter, an aeroelastic instability, involves a rotor or propeller mounted on a wing nacelle, whirling in a circle about its undeflected position.

Figure 1. Bell XV-15 tiltrotor. NASA Photo ID EC80-75, cropped. Public domain.

Aerodynamic forces acting on the blades and gyroscopic effects acting on the rotor as a whole can couple with wing modes to produce an unstable vibration which can damage or even destroy the aircraft structure. Whirl flutter claimed two L-188 Electra aircraft with all onboard, in 1959 [1] and 1960 [2], removing an entire wing in flight. It was the suspected cause of the similar loss of a Beechcraft 1900C in 1991 [3]. With their large and flexible blades, tiltrotor aircraft are particularly susceptible to whirl flutter, which generally limits their maximum cruise speed [4] as the wings must be made much thicker than would be aerodynamically optimal in order to provide sufficient stiffness. Whirl flutter remains a design consideration for turboprop aircraft [5] and is also known to affect wind turbines [6].

Existing tiltrotor designs feature gimballed hubs. The lack of bearings allows a weight saving. The flapping motion is replaced with hub tilt, which lowers Coriolis-induced loads [7] and thereby permits the use of lightweight, stiff in-plane rotor blades. These stiff blades have a first in-plane (lead-lag) frequency above one-per-rev, eliminating the possibility of ground resonance. The use of a gimbal at the hub may have a destabilising effect on the rotor’s dynamics, however. A gimballed hub is in any case an important feature of a rotor-nacelle model that is to be representative of a tiltrotor aircraft.

A great deal of literature has addressed the issue of whirl flutter in tiltrotors, using a number of approaches. One such approach is to raise the onset airspeed of whirl flutter, either to maximise the performance of an aircraft or to afford its crew a greater margin of stability at the operating airspeed. This can be achieved either by the active control of aerodynamic actuators [8, 9] or the rotor swashplate [10, 11], or by passive design measures that act against physical drivers of the whirl flutter mechanism and can thereby to an extent prevent incipient whirl flutter motion developing into larger oscillations. The first such strategy was aeroelastic tailoring [12] – the exploitation of anisotropic material properties to create dynamic couplings that are beneficial to whirl flutter stability – though rotor design modifications [13] and uncontrolled winglets [14] have also been explored. However, whatever angle the tiltrotor whirl flutter problem is approached from, at the heart of the matter lies the need for an accurate way to assess the whirl flutter characteristics of a given tiltrotor system. Despite whirl flutter being long appreciated to be a nonlinear phenomenon, from classical work such as Edenborough in the 1960s [15], through to much more recent work by Acree et al. [16], the vast majority of literature on the topic appears in almost all cases either to use analysis that is some form of linear theory, or to use time simulations that hope to observe whirl flutter behaviours directly.

Linear analyses provide localised stability predictions that do not detect changes induced by nonlinearities in the wider parameter-state space. Time simulations, while capturing the full behaviour of a nonlinear model, do not guarantee the discovery of whirl flutter solutions if they are present. A particularly important impact of nonlinearities is that they can cause the possibility of whirl flutter in parameter ranges that linear analysis declares to be stable – a phenomenon associated with what are now known as subcritical bifurcations – as shown by Breitbach [17] in another study concerning aerofoils with control surfaces. This phenomenon results in unconservative stability boundaries that can be unsafe to use [18]. This phenomenon was observed as early as 1955 by Woolston et al. [19] in work concerning the flutter of 2D aerofoils with control surfaces, though the investigation of such a phenomenon occurring in tiltrotor contexts has not been forthcoming.

Nonlinearities in a tiltrotor system may arise from a variety of sources. The drivetrain may provide both discontinuous and smooth nonlinearities [20], though other sources may include the deformability of the rotor blades and other elastic elements in the system such as elastomeric dampers [21], which linear stiffness may not adequately model. In general, the assumption of linear stiffness is only really representative of physical structures when deformations are small – a condition that may well not hold for whirl flutter oscillations – and polynomial softening and/or hardening terms may describe stiffness profiles at larger deflections more realistically [22].

Compared to the linear stability analysis methods mentioned previously, continuation and bifurcation methods (CBM) are much better suited to the stability analysis of nonlinear systems due to their ability to find solutions to these systems, and assess the stability thereof. Motivated by the phenomenon that smooth changes in system parameters can result in non-smooth changes in the stability or type of a system’s solutions, continuation is a numerical method that finds the numerical values of solutions and their respective stability as one or more system parameters are varied. Bifurcation theory then classes any changes in solution type and stability that are observed. CBM is still in the process of proliferation within the field of rotorcraft design and analysis, and its application to the tiltrotor whirl flutter problem appears to be a novel idea.

The authors investigated a basic non-gimballed rotor-nacelle system in 2017 [23] and 2018 [18], where a smooth nonlinearity was introduced into the shaft yaw stiffness and the impacts of the nonlinearity were summarised and discussed using stability boundaries. A brief investigation into the impacts of such a smooth nonlinearity on a gimballed hub model was conducted in 2019 [24], though the authors have since discovered new results that constitute significant changes to the stability boundary and therefore warrant attention. Furthermore, some complex behaviours were discovered in supposedly stable parameter regions that are analysed here. The aim of this study is to explore and understand the impact of the aforementioned smooth nonlinearity on a gimballed rotor-nacelle model. Through this approach, the importance of using more representative nonlinear structural stiffness models in analysing whirl flutter stability of tiltrotor aircraft will be shown.

In this paper, the rotor-nacelle model, including details of the nonlinear adaptation, is presented in Section 2. The stability analysis methods used in this paper are described in Section 3. These are applied in Section 4: linear stability analysis is used to establish the baseline whirl flutter stability characteristics of the system, followed by continuation and bifurcation analysis to uncover the effects of the nonlinearity’s presence. Stability boundaries are used to summarise these effects. Further analysis of periodic and quasi-periodic behaviours, via Poincaré sections and power spectra, is also presented in this section.

2.0 Whirl flutter model

2.1 9-DoF whirl flutter model

An 18-state, 9-DoF model originally presented by Johnson [25] 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 [26]. 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. 25.

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 [25], 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. 24.

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. 25 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.0 Stability analysis methods

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 [28], 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 [29].

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 [30], which uses an implementation of AUTO-07P [31]. 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. 32.

4.0 Results and discussion

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 [25]. 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. 25, 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 [32]. 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.

5.0 Conclusions

This article has demonstrated the use of continuation and bifurcation methods to provide nonlinear dynamic analysis of whirl flutter. Polynomial stiffness profiles are frequently more representative of structures over large deflections than the linear profiles that are sometimes used in models in tiltrotor whirl flutter research.

To investigate the impacts of such nonlinear stiffness on tiltrotor whirl flutter predictions, a well-established linear rotor-nacelle system model was implemented and adapted to feature the nonlinear stiffness profile in the rotor blades’ flapping stiffness. The nonlinear stiffness profile used was a cubic polynomial function with a positive cubic coefficient to provide hardening characteristics. Appropriate stability analysis methods were described and employed for both the linear and nonlinear models, with continuation and bifurcation methods (CBM) being used in the latter case. Eigenvalue analysis was used to establish a baseline stability boundary for the model, representing both the linear model version’s stability characteristics, and also what the linear analysis “sees” in the nonlinear version of the model. Bifurcation diagrams were generated for a number of wing torsion stiffness cases for the model.

The results showed that whirl flutter was possible in a parametric region where linear analysis predicts local stability. That is, stable whirl flutter LCOs were found to exist in this region, giving the possibility of the system joining them following a sufficient perturbation, such as a gust or loads induced by manoeuvring. This phenomenon of stable whirl flutter LCOs existing in supposedly stable regions of the stability boundary was termed “overhang” on account of its appearance on bifurcation diagrams.

The stability boundary was redrawn to take account of the overhang phenomenon, which appended an extra “unsafe” region to the existing unstable region. Important of note is that the prediction of the overhanging LCO to be stable is not an indication that encountering it is safe, but rather that this LCO can attract the system to it if it is perturbed appropriately. Allowing a real-world tiltrotor system to encounter it should not be thought of as viable. Upon entering the LCO, the large oscillation amplitudes would likely cause a rapid degradation in the structural properties of the aircraft, most likely leading to structural failure. Even if this degradation were not to occur immediately, the oscillations would present a fatigue hazard to aircraft nacelle mounts.

Among the overhanging whirl flutter branches, the torus bifurcations and the quasi-periodic behaviour that they entail are likely a still greater threat than the LCOs. The larger amplitude and broadband-frequency motions that they involve are likely to accelerate the degradation of the tiltrotor system’s structure, leading either to accelerated fatigue in the best case or expedited catastrophic structural failure in the worst.

While several of the predicted motions are likely too large to exist for a full cycle, let alone form a persisting motion, the whirl flutter insights gained – particularly those concerning the type of solution (LCO/quasi-periodic) are still of use in systems where smaller amplitude solutions occur. Additionally, the consideration of structural strengths in modelling would allow motions causing structural failure to be identified with confidence, rather than by assumption.

Future improvements of the model that could be undertaken include refinements to the aerodynamic model, and the further development of the rotor blade dynamics to include coupled mode shapes and a torsional degree of freedom.