Application of continuation and bifurcation methods to aeroelastic blade stability

Paper number : 132

Application of Continuation and Bifurcation Methods to Aeroelastic Rotor

Blade Stability

Djamel Rezgui

AgustaWestland/Bristol University Technology Centre

Dept. of Aerospace Engineering University of Bristol, UK djamel.rezgui@bris.ac.uk

Mark H. Lowenberg

Dept. of Aerospace Engineering University of Bristol, UK m.lowenberg@bris.ac.uk

Mark Jones

AgustaWestland Ltd, Yeovil, Somerset, UK mark.jones@agustawestland.com

Claudio Monteggia

AgustaWestland SPA, Cascina Costa, Italy claudio.monteggia@agustawestland.com

Abstract

The present paper illustrates how continuation and bifurcation methods are applied to investigate the aeroelastic stability of helicopter blades over a range of flight and control parameters. The paper first compares the bifurcation results for flap-lag instability to those available in the literature, which are obtained using conventional methods. In this case, the continuation analysis is applied to a simple rotor system incorporating a rigid blade with flap and lag degrees of freedom with root springs to model the stiffness of hinges. The continuation analysis is then carried out to investigate the aeroelastic stability of rotor blades, where modal representation of the blade flexibility is used. The analysis is performed in hover and forward flight, within and beyond operating flight conditions. The significance of the typical bifurcations that the blade behaviour undergoes is discussed.

N

′ _{Differentiation with azimuth angle.}

(˙) Differentiation with time.

Alat,Blong Lateral and longitudinal cyclic pitch angles.

CT Rotor thrust coefficient.

CL, CD, CM Elemental lift, drag and pitching moment coefficients.

Cd0 Blade profile drag coefficient.

L, D Elemental aerodynamic flapwise, lagwise force com-ponents, positive upward and forward respectively.

L0, D0 Elemental aerodynamic flapwise and lagwise force components computed when the blade’s modes are evaluated, positive upward and backward respectively.

M Elemental aerodynamic pitching moment about the section shear centre, positive nose up.

M0 Elemental aerodynamic pitching moment computed when the blade’s modes are evaluated, positive nose up.

Ma Elemental Mach number.

NB Number of blades.

Qqi Generalised force of mode i.

R Rotor radius.

T Rotor oscillation period.

Ti, Tav Rotor instantaneous and average thrust.

Treq Required rotor thrust.

UP, UT Perpendicular and tangential components of the ele-mental flow velocity.

Ω Rotor speed.

Θ Blade pitch due to twist deformation.

α Elemental angle of attack.

αs Rotor shaft angle. ¯

R Elastic coupling parameter. ¯

f Helicopter flat plate drag area.

β Flap angle, positive upward.

βa1,req,βb1,req Required lateral and longitudinal flapping angles.

βpc Precone angle.

δr Width of the blade element.

δ Parameter vector describing the rotor and flow proper-ties.

η negative real part of the eigenvalue associated with the lag degree-of-freedom.

γ Lock number.

λ Inflow ratio.

µ Advance ratio.

ν0,νs,νc Average, lateral and longitudinal induced velocity components.

νi Induced velocity at position r and azimuthψ.

ωβ,ωζ Non-rotating flap and lead-lag frequencies.

ωqi Modal frequency of mode i.

φ Inflow angle.

ψ Azimuth angle.

σ Rotor solidity.

θcol Collective pitch angle.

θpre.twist Blade local built-in twist.

ζ Lag angle, positive forward.

c Blade elemental chord.

p Dimensionless rotating flapping frequency atθ= 0,

p =

q 1 +ω2_β.

qi Generalised displacement of mode i.

r Elemental radial position.

t Time.

u, y Harmonic oscillator states, u = sin(Ω t) and y =

cos(Ωt).

wi,vi,ti Flap, lag, and torsion mode shapes of mode i.

x State vector.

Ii Modal mass of mode i.

I

The stability of a helicopter rotor blade is a complex dynamical problem, involving aerodynamic, structural, material and geometric nonlinearities. Furthermore, the trend towards higher performance gains and lower vibra-tion and noise levels has led to the development of more complex rotor systems, which incorporate novel designs features, utilizing for example: composite materials, ac-tive lag dampers and acac-tive trailing edge flaps. These features tend to increase the levels of nonlinearity in the rotor system, which means that a proper nonlinear analy-sis of the blade dynamics is required. Until this time, dif-ferent mathematical techniques have been used to study the aeromechanical and aeroelastic blade stability, at dif-ferent flight regimes. These techniques include time his-tory simulation (time integration techniques) [1], para-metric resonance analyses [2, 3, 4, 5], perturbation meth-ods [6, 7] and Floquet analysis [7, 8, 9]. In fact, the he-licopter blade aeromechanical and aeroelastic stability is well understood in both the academic and industrial sec-tors, including cases of very high tip speed ratios. Exten-sive reviews in the field are those by Friedmann [10, 11], Bielawa [12] and Chopra [13].

However, many of the above stability methods de-pend on assumptions that are questionable for newer ro-tor configurations, or may not provide the complete sta-bility picture. For example, the methods can predict the local stability of the blade, but the regions of attraction in that case are not defined. In other words, the blades might be stable for small disturbances but not necessar-ily for large ones; and hence the outcome in the event of a large disturbance inducing instability is not indicated.

In recent years, the stability of rotors in autorotation was investigated by Rezgui et al. and Lowenberg et al. using nonlinear dynamics theory implemented numeri-cally in the form of continuation and bifurcation meth-ods [14, 15, 16, 17]. The studies showed that these tech-niques are crucial in the identification of the instability scenarios of rotors in autorotation including but not ex-clusive to blade sailing, high speed instability, control input-induced instability and lightly loaded rotor insta-bility. The analysis showed the presence of different bi-furcation points which establish the coexistence of stable and unstable periodic solutions over a range of param-eters. The analysis was complemented by wind tunnel testing, the results of which exhibited substantial agree-ment with those obtained from bifurcation analysis, in-cluding an unexpected asymmetric form of stable autoro-tation.

In the present paper, the continuation and bifurcation methods are applied to investigate the aeroelastic stabil-ity of a helicopter blade over a range of flight and control parameters. The analysis first investigates the stability of a rigid blade model and compares the bifurcation results to those available in the literature, which are obtained using conventional methods. In this case the blades are modelled to incorporate rigid lag and flap degrees of free-dom with root springs to represent the stiffness of hinges. The continuation analysis was then carried out to investi-gate the aeroelastic stability of rotor blades, where modal representation of the blade flexibility is used.

C

B

M

B

S

The basic idea of the continuation techniques is the cal-culation 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 so-lutions construct a number of equilibrium branches that could be either stable or unstable. To determine the sta-bility, either an Eigen or Floquet analysis is carried out at each computed solution, depending on the nature of this solution. For instance, in hover the blade behaviour can be said 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, hence Flo-quet theory is used to determine the stability.

Bifurcation is the qualitative change in the system be-haviour 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 multi-ple solutions for the same set of system parameters. The identification of these different solution branches helps to uncover the global dynamics of the system. Particu-lar 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 bifurca-tion points. Then, follow the emerging branches to con-struct a more complete picture of the system dynamics (bifurcation diagram). The different types of bifurcations that can occur in equilibria or periodic orbits are not dis-cussed in this paper. However, the reader is referred to general texs such as references [18, 19] for more back-ground on the subject.

The continuation algorithm used in this analysis is implemented in the continuation and bifurcation soft-ware AUTO [20]. Besides many other types of equa-tions, AUTO can perform extensive bifurcation analysis of ordinary differential equations (ODEs) of the form:

˙

x(t) = f (x(t),δ), f, x ∈ℜn (1)

subject to initial conditions, boundary conditions, and in-tegral constraints. Here x is the state vector and δ de-notes one or more parameters. Equation (1) is written in the generic (nonlinear) space form, where the state-derivatives are functions of the states and some param-eters. Therefore, as far as AUTO is concerned, the dy-namical model must be constructed in a format which al-lows passing the state derivatives to AUTO and receiving states and parameters from it, regardless of the environ-ment where the system equations are coded. This soft-ware attribute provides a stability tool which not only is independent from the model but also able to couple with a wide range of modelling platforms, which can be par-ticularly beneficial in the helicopter industry.

D

N

R

B

M

In this paper, the continuation and bifurcation analysis is applied to two different rotor blade models. Both are modeled in the dynamical system form of Equation (1) where the parameter vectorδcan hold values for blade

control angles, helicopter velocities, blade properties, etc. Furthermore, the ODEs in Equation (1) are au-tonomous in that the independent variable t does not ap-pear explicitly in the equations. However, in cases where the system is periodically forced, as in the dynamics of rotor blades in forward flight, the independent variable t (or the azimuth angleψ) has to be converted to a state variable. Of course, one can always designate the time t or azimuth angleψas additional states, in order to trans-form the system to an autonomous one. This can either be achieved by ˙ψ=Ω, where ψ is a state or by ˙t = 1, where t is a state and the azimuth angle can be calculated

asψ=Ωt sinceΩcan be assumed constant. However,

if the above method is used (either cases), the new time state will monotonically increase and hence will not de-scribe an oscillatory behaviour. Therefore, to solve this problem, a harmonic oscillator model can be used to re-alise the periodicity of all states. The harmonic oscillator equations are: ˙ u = u + Ωy − u (u2_{+ y}2₎ ˙ y = −Ωu + y − y (u2_{+ y}2₎ (2) or u′ = u + y − u (u2_{+ y}2₎ y′ = u + y − y (u2_{+ y}2₎ (3)

where the u = sin(Ωt) and y = cos(Ωt) are solutions

to Equation (2). The terms u and y can now be used to replace any sin(Ω t), sin(ψ), cos(Ω t) or cos(ψ)

in the blade forcing equations as appropriate. Alter-natively, the azimuth angle can be calculated using the quadrant-arctangent function (atan2 function in MAT-LAB or FORTRAN).

ψ= atan2(u, y)

Rigid Blade Model with Flexible Hinges

In this section, continuation and bifurcation methods are applied to investigate the flap-lag stability of a rigid he-licopter blade with springs attached between the blade root and rotor shaft, in the flap and lag degrees of free-dom. This spring model is sufficient to represent an elas-tic blade and hub with hinge offsets. The model used here is published by Peters [7].

Equations of Motion

The equations of motion of the rotor blade in the flap and lag degrees of freedom are given by:

β′′_{+ sinβcosβ(1 +ζ}′₎2_{+ (P− 1)(}β₋β

pc)

+ Zζ=_Ω12_I

RR

0 Fβrdr

cos2βζ′′− 2 sinβcosβ(1 +ζ′)β′+ Wζ + Z(β−βpc) =cos(_Ω2β)_I

RR

0 Fζrdr

where P, Z and W are stiffness parameters which depend on the non-rotating flap and lead-lag frequenciesω_βand ω_ζ. F_βand F_ζ are the force components per unit length perpendicular and parallel, respectively, to the direction of rotation. The expressions for the stiffness parameters and the force components are given in reference [7].

Rotor Model in the State Space Form

The main model equations, which have been generated in Equations (3) and (4) can be rearranged in terms of time derivatives of the state vector x to produce the required state-space form for the stability analysis:

x′=                u y β′ β′′ ζ′ ζ′′                = f (x,δ) where x =                u y β β′ ζ ζ′                (5)

whereδ={µ,λ,θcol, Alat, Blong,ωβ,ωζ, . . .}T.

Aeroelastic Rotor Blade Model

In this section, modal representation is used for the struc-tural blade dynamics and hence each blade is represented by a number of general modes (eight in this analysis). The mode shapes, frequencies, and modal masses are computed a priori either in vacuum or at a given hover case condition. Although the model can include the dy-namical equations for all blades of the rotor, it is ade-quate here to use a single bladed rotor model for the sta-bility analysis.

The forced response equation approach is used to de-scribe the aeroelastic dynamics of the blade. The equa-tions used were developed in reference [21]. The formu-lation of this equation was done in a manner such that the orthogonality of the modes led to an equation in which the modal response depends only on the following forc-ing components:

1. the aerodynamic forcing;

2. time dependent terms such as the Coriolis force; 3. blade pitch dependent terms;

4. nonlinear terms not included in the formulation of the modal equation.

For the aerodynamic forcing, a blade element tech-nique is used for calculating the forces and moments act-ing on the blades. This technique is based on dividact-ing each blade into a number of elements, then the aerody-namic forces are computed for each element, consider-ing them as quasi-two-dimensional aerofoils, with as-sociated lift, drag and pitching moment characteristics.

To achieve this, elemental velocity components need to be determined. Unsteady aerodynamics, blade/blade and blade/airframe interactions are not considered here. To model the inflow, a three state dynamic induced veloc-ity model is used. The code of this model is written in MATLAB.

Modal Response Equation

The forced response equation for each mode is: ¨ qi+ω2qi qi= Qqi Ii i = 1, 2, 3, ... (6) hence ¨ qi=−ω2qi qi + Qqi Ii i = 1, 2, 3, ... (7)

whereωqi, Ii, qi and Qqi are the modal frequency, modal mass, generalised displacement and generalised force terms for each mode i. The differentiation ( ¨ ) is done with respect to time t. It can be noticed from equation (6) that the elastic damping is assumed negligible and hence is not accounted for. Each mode shape consists of flap, lag and twist mode shape components, wi, vi and ti

re-spectively. The generalised force Qqi consists of the four forcing terms described above and can be written as fol-lows:

Qqi = Qaero,qi+ Qtdep,qi+ Qpper,qi+ Qnlin,qi (8)

where the expressions for the forcing terms can be found in [21].

Local Flow Velocities

The blade deformation due to bending and twisting af-fects the local flow velocities and also angles of attack. In general, the components of the resultant flow velocity at each blade element arise from five sources, namely: blade rotation, free stream due to the helicopter move-ment, rotor induced velocity, rates of blade bending and rates of change of the predeformed blade coordinates. Furthermore, the resultant flow velocity is convention-ally resolved into components tangential (UT) and

nor-mal (UP) to the local axes of the blade. The full

expres-sions of these components depend on many variables, in-cluding the positions of the blade elements before and after deformation, and can be very long. Hence, they are not presented here but can be found in [21].

Blade Aerodynamics

The aerodynamic forcing Qaero,qi in Equation (8) is the most complex term to evaluate in comparison to the rest

of the forcing terms. This is due to the complicated na-ture of the flow field around the aerofoil section. The general form of Qaero,qi is given by the integral:

Qaero,qi = RR 0 h d(L−L0) dr wi+ d(D+D0) dr vi +d(M−M0) dr ti i dr (9)

where the three terms on the right hand side represent the force distribution in the flapwise and lagwise direc-tion and the pitching moment distribudirec-tion respectively.

To evaluate the above aerodynamic loads, the ele-mental lift and drag forces and pitching moment need first to be calculated. In this analysis, the blade ele-ment method is adopted an hence two-dimensional quasi-steady flow is assumed. The expression for the elemental lift, drag and pitching moment are given as:

δL = 1 2ρU 2_cδrC L(α, Ma, r) (10) δD =1 2ρU 2_cδrC D(α, Ma, r) (11) δM =1 2ρU 2_c2_δrC M(α, Ma, r) (12) whereρ, c, U = q

U_P2+ U_T2,δr, CL, CDand CMare the

lo-cal air density, the blade elemental chord, the elemental resultant flow velocity, the width of the blade element, and the lift, drag and pitching moment coefficients at each element respectively. α and Ma are the elemen-tal angle of attack and Mach number respectively. The aerodynamic force and moment coefficients for each in-dividual blade element are calculated numerically using nonlinear look-up tables. Experimental data for chosen aerofoil sections are used. These tables provide data for a 360◦range of angle of attack and local Mach numbers of up to 0.85. The local angle of attack is given by:

α=θ+φ+θquasi (13)

whereθis the elemental pitch angle, θquasi= _2Uc˙θ is the

aerodynamic quasi-steady effect on pitch angle andφis the elemental inflow angle, which can be calculated as follows:

φ= atan2(UP,UT)

The local pitch angle θis a combination of all the local pitch angle contributions, i.e.:

θ=θcol− Alatcos(ψ)− Blongsin(ψ) +θpre.twist+Θ

whereθcol, Alat,Blong,θpre.twist and Θare the blade

col-lective pitch, the lateral cyclic, the longitudinal cyclic, the blade local built-in twist and the blade pitch due to twist deformation.

Finally, the local forces acting on a blade element in the blade coordinate system can therefore be determined from the elemental lift and drag forces.

L =δL sinφ−δD cosφ

D =−δL cosφ−δD sinφ (14)

whereas, the local elemental aerodynamic pitching mo-ment is simply defined as:

M =δM (15)

The rotor instantaneous thrust Ti can be evaluated

simply by summing all elemental vertical forces L and multiplying it by the number of blades NB. i.e.

Ti= NB

N

∑

elem=1

L (16)

This value of thrust is only used to estimate the induced velocity within the rotor model. The correct thrust value, which is used for performance and trimming procedures, has to be averaged out across one rotor revolution.

Tav= 1 2π Z 2π 0 Tidψ (17)

Induced Velocity Model

The inflow is captured via a 3-state Pitt-Peters dynamic wake model [22, 23, 24]. This model permits the vari-ations of the induced velocity in both the radial and az-imuthal position. Furthermore, it allows the lag dynam-ics associated with moving a volume of air to be mod-elled. The inflow model is given for three states as fol-lows:

νi(r,ψ) =ν0+

r

R(νssinψ+νccosψ) (18)

whereνi is the induced velocity at an element of radius

r and azimuth positionψ. The induced velocity

compo-nentsν0,νsandνcare given in the wind axes by:

[τ]   νν˙˙0s ˙ νc   w =−   ν_ν0_s νc   w + [L]   LTaeroaero Maero   w (19)

Taero, Laero and Maero are the thrust, the aerodynamic

rolling and pitching moments respectively in the wind axis and expressions for the matrices [τ] and [L] can be

found in [22, 23, 24].

Rotor Model in the State Space Form

The main model equations, which have been generated in Equations (2), (7) and (19), can be rearranged in terms

of time derivatives of the state vector x to produce the required state space form:

˙ x =                                  ˙ u ˙ y ˙ ν0 ˙ νs ˙ νc ˙ q1 ¨ q1 .. . ˙ q8 ¨ q8                                  = f (x,δ) where x =                                  u y ν0 νs νc q1 ˙ q1 .. . q8 ˙ q8                                  (20)

C

B

R

Rigid Blade Model with Flexible Hinges

In order to find the steady state or periodic solutions for the rigid blade model as given in Equation (5), the inflow ratioλand the control anglesθcol, Alat and Blongneed to

be known for a given set of µ, ω_β, etc. Of course, the analysis can be carried out at a prescribed values ofλ, θcol, Alat and Blong. However, the conventional way is

to computeλbased on a momentum theory assumption, and find the control angles which satisfy the intended trim requirements.

Unlike many analyses, where rotor trimming is car-ried out first (using the harmonic balance method for example) prior to any stability calculations, AUTO can compute solutions and their stability simultaneously for given boundary or integral conditions. This means that the continuation and trimming procedure can be done in parallel. In fact, the eigenvalues and Floquet multipli-ers are computed at negligible extra cost to the continua-tion analysis. To findλ,θcol, Alat and Blongfour integral

boundary conditions are constructed as follows:

1. λis obtained by equating the thrust calculated by integrating the vertical force along the blade to the thrust from simple momentum theory:

2νpµ2₊λ2₋ σa 2πγcosβ Z 2π 0 1 Ω2_I Z R 0 F_βrdrdψ= 0 (21) where the induced velocityνis given by:

ν=λ− µαs

and αs is the rotor shaft angle, which can be

ob-tained from the propulsive time condition αs=

µ2f¯

2CT

where ¯f is the helicopter flat plate drag area. ¯f can

be set to zero if propulsive trim is not required.

2. θcol is obtained by equating the calculated thrust

obtained by integrating the vertical force along the blade to the required thrust:

CT,req− σ a 2πγcosβ Z 2π 0 1 Ω2_I Z R 0 Fβrdrdψ = 0 (22) 3 & 4 To find Alat and Blong, moment trim condition is

assumed at the rotor hub. This condition is ob-tained by suppressing the first harmonic compo-nents of the flapping angle using the cyclic pitch angles. This is implemented as follows:

Z 2π 0 βsinψdψ= Z 2π 0 βu dψ= 0 (23) Z 2π 0 β cosψdψ= Z 2π 0 β y dψ= 0 (24)

It should be noted that high harmonic components of the blade flapping angle are not suppressed by Equations (23) and (24). 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −3 −2 −1 0 1 2 3x 10 −3 Advance ratio µ η f−=0 f−=0.01 STABLE UNSTABLE

Figure 1: Lag damping variation with advance ratio. p = 1.15,

ωζ= 1.4,R = 0¯ ,γ= 5,σ= 0.05,Cd0= 0.01,CT/σ= 0.2.

The rigid blade model was coded in MATLAB, and AUTO was supplied with the state-space equations of the model and the four integral conditions, i.e. equa-tions (5), (21), (22), (23) and (24). The continuation was performed at two cases of trim condition. The first is a propulsive trim at ¯f = 0.01 and the second is a

mo-ment trimmed condition in which ¯f = 0. In order to

compare the continuation results with those published by Peters [7], the damping parameterηassociated with the lag degree of freedom is extracted from the computed Floquet multipliersϑ. To findηthe logarithm of the Flo-quet multiplier associated with the lag mode is evaluated to find the equivalent eigenvalue. Then ηis simply the negative real part of this eigenvalue (η=−ℜ(lnϑ)).

Figure 1 shows the variation of lag damping with ad-vance ratio for ¯f = 0 and ¯f = 0.01. Like the predictions

unstable at certain ranges of advance ratios. The damp-ing graph for the case ¯f = 0.01 is identical to the graph

predicted by Peters. However, for the case ¯f = 0,

al-though ηbecomes negative at the same µ predicted by Peters (≈ 0.16), the value at which η becomes positive again is higher when using the continuation method (≈

0.73 compared to ≈ 0.49). This discrepancy is thought to

be due to the different trimming procedures. Peters used the harmonic balance method where second and higher harmonics of the blade flapping were assumed negligi-ble, whereas no such assumptions were used in the con-tinuation analysis. In fact, although the first flap harmon-ics were suppressed to ensure moment trim requirement, it can be shown that the second flap harmonics are large for µ> 0.4, when the trimming is achieved using the

con-tinuation approach. Therefore, the values of the control angles evaluated by the two approaches can be different.

0 0.2 0.4 0.6 0.8 1 5 10 15 20 TR TR Advance ratio µ

Peak flapping angle (deg) 50 0.5 1 1.5

10 15 20 25 30 35 TR TR PD PD Advance ratio

Peak flapping angle (deg)

0 0.2 0.4 0.6 0.8 1 −1 0 1 2 3 4 5 TR TR Advance ratio µ

Peak Lag angle (deg)

0 0.5 1 1.5 −5 0 5 10 15 20 TR TR PD PD Advance ratio µ

Peak lag angle (deg)

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 TR TR Advance Ratio µ Inflow ratio λ (a) ¯f = 0.01 0 0.5 1 1.5 0 0.02 0.04 0.06 0.08 TR TR PD _PD Advance ratio µ Inflow ratio λ (b) ¯f = 0

Figure 2: Bifurcation diagrams for the cases f = 0¯ .01 and

¯

f = 0. Continuation parameter = µ. Solid line: stable

peri-odic branch, dashed line: unstable periperi-odic branch. For the flap and lag angles, only the peak values of the oscillations in one cycle are plotted. TR denotes torus bifurcations.

Figure 2 depicts the bifurcation diagrams for the cases of ¯f = 0.01 and ¯f = 0 respectively. The variation

of inflow ratioλcomputed by AUTO is also plotted for completeness. The bifurcation diagrams show that the change in stability of the periodic branches, as previously predicted by η changing sign in Figure 1, corresponds to the presence of torus (Neimark-Sacker) bifurcation points. This type of bifurcation occurs when a pair of complex conjugates of the Floquet multipliers cross the unit circle (in the real and imaginary plane) at a non zero

real part. Furthermore, a new secondary quasi-periodic branch (branch of invariant tori) emerge from the bifur-cation points. To investigate this further, the rotor model was simulated in time from the unstable branch at µ = 0.5

for the case of ¯f = 0. The appropriate values forλ,θcol,

Alatand Blongwere used. Figure 3 illustrates that the

flap-ping angle oscillation starts from the unstable periodic solution (limit cycle) and then diverges to a stable quasi-periodic solution. The unstable quasi-periodic solution has a time period of 1 rev (2πrad), while it is not possible to identify an oscillation period for the stable quasi-periodic behaviour of the blade flapping motion. To further ver-ify that the solution is quasi-periodic, a segment of the time response was plotted in both (ψ,β, ˙β) and (ψ,ζ, ˙ζ) cylindrical coordinate systems. It can be shown that so-lution trajectories trace the surface of a torus but do not connect. 0 200 400 600 800 1000 1200 1400 1600 −25 −20 −15 −10 −5 0 5 10 15 20 25 30 Azimuth (rev)

Flapping angle (deg)

0 1 2 3 2 4 6 8 10 Azimuth (rev)

Flapping angle (deg)

1560 1570 1580 1590 −20 −10 0 10 20 30 Azimuth (rev)

Flapping angle (deg)

Figure 3: Time simulation from unstable branch atµ = 0.5and

¯

f = 0.p = 1.15,ω_ζ= 1.4,R = 0¯ ,γ= 5,σ= 0.05,Cd0= 0.01,

CT/σ= 0.2.

(a) Flap DOF (b) Lag DOF

Figure 4: Segment of the stable quasi-periodic solution traces a surface of a torus in a cylindrical coordinate frame of axes.

µ = 0.5, f = 0¯ , p = 1.15,ω_ζ= 1.4,R = 0¯ ,γ= 5,σ= 0.05,

Cd0= 0.01,CT/σ= 0.2.

Two more characteristics of the unstable periodic branch can be observed from Figure 3. First, it can be seen that the second harmonics are evident in the flap-ping motion and are not suppressed. Second, the time taken for the flapping motion to clearly diverge from the

unstable periodic orbit toward the quasi-periodic solution is very high (about 400 revolutions), which may lead the solutions to appear to be stable. This means that a good level of accuracy in the stability method is necessary to predict the correct stability for this system at this param-eter combination. Figure 2-(b) also illustrates the pres-ence of period doubling bifurcations at high advance ra-tios. The first bifurcation occurs at µ≈ 0.9 where the

main stable periodic branch becomes unstable, whereas the second bifurcation is at an even higher value of µ. Pe-riod doubling bifurcation occurs when one Floquet multi-plier equals−1 with zero imaginary part. After the

bifur-cation point, the main periodic branch changes stability and a new secondary periodic branch emerges. The peri-odic solutions of this secondary branch have double the period of the solutions of the main branch. To illustrate this point, the rotor model was simulated in time from the unstable branch at µ = 1 for the case of ¯f = 0. The

appro-priate values forλ, θcol, Alat and Blongwere used.

Fig-ure 5 illustrates that the flapping angle oscillation starts from the unstable periodic solution (limit cycle) and then diverges to another stable periodic solution. The unstable main periodic solution has a time period of 1 rev, while the other secondary periodic solution has a time period of 2 rev (3 cycles between 1584 and 1590 rev).

0 200 400 600 800 1000 1200 1400 1600 −80 −60 −40 −20 0 20 40 Azimuth (rev)

Flapping angle (deg)

0 1 2 3 −10 0 10 20 30 Azimuth (rev)

Flapping angle (deg)

1584 1586 1588 1590 −40 −20 0 20 40 Azimuth (rev)

Flapping angle (deg)

Figure 5: Time simulation from unstable branch atµ = 1.0and

¯

f = 0.p = 1.15,ω_ζ= 1.4,R = 0¯ ,γ= 5,σ= 0.05,Cd0= 0.01,

CT/σ= 0.2.

Aeroelastic Rotor Blade Model

Similar to the rigid blade analysis, the rotor trim condi-tions were also embedded within the continuation proce-dure. Since a dynamic inflow model is used in this case, only three integral conditions were needed to be speci-fied. These can be described as follows:

Tav= 1 2π Z 2π 0 Ti dψ= Treq (25) 1 π Z 2π 0 β sin(ψ) dψ=1 π Z 2π 0 β u dψ=βa1,req (26) 1 π Z 2π 0 β cos(ψ) dψ=1 π Z 2π 0 β y dψ=βb1,req (27)

The rotor model was configured to have realistic blade and flight characteristics. Figure 6 illustrates the continuation results when the advance ratio was used as the continuation parameter. Only the peak values of os-cillatory modal displacements of the first two modes are plotted. The bifurcation diagrams show that when µ is between approximately 1.03 and 1.19 the periodic branch is unstable, due to the presence of two torus bifurca-tions. To further scrutinise the results, the variation of the damping of the modes was investigated from the com-puted Floquet multipliers. It was found that the unstable periodic branch corresponds to the range of advance ra-tios when the damping of the first lag mode (mode 1) is negative. Hence, this indicates that this instability may be of the flap-lag type similar to the one obtained in the previous section. 0 0.5 1 1.5 −50 0 50 100 Advance ratio µ q1 TR TR 0 0.5 1 1.5 0 100 200 300 400 500 Advance ratio µ q2 TR TR

Figure 6: Bifurcation Diagram for forward flight case.

Continu-ation parameter =µ. Solid line: stable periodic branch, dashed

line: unstable periodic branch. Only the peak values of the oscillations in one cycle are plotted. TR denotes torus bifurca-tions.

It should noted that although the model might not be suited for flight conditions of very high advance ratios (µ> 0.5 for example), it is customary when performing

continuation analysis to extend the continuation param-eter beyond the physical range. The reason for this is to search for any bifurcation points that might lead to new solution branches, which return back to the physi-cal range. For the case investigated in Figure 6, the only secondary branch found was the quasi-periodic branch, for 1.03 ≤ µ ≤ 1.19 emerging from the torus bifurcations

points.

Continuation runs were also carried out in the hover case over a range of collective pitch angles. The resulting bifurcation diagrams are depicted in Figure 7. In stable hover condition, the blade is in steady state situation and hence its behaviour is not periodic. For this reason there is no need to use neither the harmonic oscillator method, since there is no periodic forcing to the system, nor to trim the rotor to certain flapping angle components. Fur-thermore, since the collective pitch is used as the contin-uation parameter, which will directly affect rotor thrust, the integral thrust condition is not necessary for this case.

0 5 10 15 20 25 30 35 40 −1200 −1000 −800 −600 −400 −200 0 200

Collective pitch angle (deg) q1

26 27 28 −700

−650 −600

Collective pitch angle (deg)

q1 HB LP 20 21 22 23 −350 −300 −250 −200 −150 HB

Collective pitch angle (deg)

q1 HB HB HB HB HB LP LP LP TR LP LP TR TR 0 5 10 15 20 25 30 35 40 −100 0 100 200 300 400 500 600 700

Collective pitch angle (deg) q2 20 21 22 23 50 60 70 80 90 100 110 120 HB

Collective pitch angle (deg)

q2 HB HB HB HB HB LP LP LP TR LP TR LP TR 27 27.1 27.2 27.3 200 205 210 215

Collective pitch angle (deg)

q2

HB LP

Figure 7: Bifurcation Diagram for hover case. Continuation

pa-rameter =θcol. Black solid line: stable periodic branch, black

dashed line: unstable periodic branch, gray solid line: sta-ble equilibrium branch, dotted gray line: unstasta-ble equilibrium branch. Only the peak values of the periodic oscillations in one cycle are plotted. TR, HP and LP denotes torus, Hopf and limit point (fold) bifurcations respectively.

Asθcol is increased from 0◦, there is only one stable

equilibrium (hover) branch, which becomes unstable be-tween 19.8◦_{and 23}◦_{. This equilibrium branch also}

dista-bilises atθcol higher than 23◦. In the first instability case,

the main branch changes stability at two Hopf bifurcation points located at 19.8◦ _{and 23}◦_{. This type of}

bifurca-tion occurs when a complex conjugate pair of eigenval-ues cross the imaginary axis, of the real and imaginary plane, with non zero imaginary parts. Hopf bifurcation is also associated with the birth of secondary periodic solutions. In Figure 7, both Hopf points are subcritical leading to the birth of unstable periodic solutions, which co-exist with sections of the stable equilibrium branch. Subcritical Hopf is a hard bifurcation which will lead to the solutions to jump to another attractor just after the bi-furcation point. This jump can be very hazardous. The unstable periodic branch emerging from the Hopf point atθcol = 19.8◦ extends very slightly as θcol is reduced

and then it folds back at a limit point (fold bifurcation) atθcol = 19.6◦. This bifurcation point occurs when one

Floquet multiplier crosses the unit circle at a value of 1. The solution then becomes stables and extends as θcol

is increased. Furthermore, this stable periodic branch changes stability a few times in the vicinity ofθcol= 23◦

due to the existence of three other fold bifurcation points, until it merges with the main equilibrium branch at the Hopf point located at θcol = 23◦. It can be shown that

if the blade is disturbed from the unstable equilibrium branch, it will always get attracted to the stable periodic branch. In other words, if the pitch angle is between 19.8◦and 23◦, the blade will eventually behave in a peri-odic manner.

The Hopf point at θcol = 19.8◦ is found to be

asso-ciated with the first flap mode becoming unstable, which stabilises again atθcol= 23◦. It was also found that other

two flap modes lose and gain stability within the same range of pitch angles. This is depicted by other four Hopf bifurcation points within the unstable equilibrium branch. The unstable periodic branch which emerge from two of these points is plotted in Figure 7. This branch also experience torus bifurcations.

Figure 7 also illustrates the existence of a subcriti-cal Hopf and limit point bifurcation, which fold the sec-ondary unstable periodic branch into a stable periodic one. If the pich angle is increased beyond 27.23◦_{, the}

solutions will jump to the stable periodic branch. It can be also seen that in the range θcol = 27.07◦ to 27.23◦

the blade can behave either in a stable equilibrium or in a stable periodic manner, depending on the perturbation levels subjected to the blade. Although this pitch range is small and non-operational, it provides a good example on how important the construction of the bifurcation di-agram is in predicting the global nonlinear dynamics of the blade behaviour.

C

This paper presented an implementation of continuation and bifurcation methods in studying the aeroelastic sta-bility of helicopter blades over a range of advance ratios and blade pitch angles. The set up for the continuation analysis and the general blade modelling layout were de-scribed for two different nonlinear blade models. These aeroelastic models were constructed in a generic state space form, where the harmonic oscillator approach was used to transform the periodically forced equations into autonomous ones. Furthermore, the paper described how the trimming procedure is implemented within the tinuation analysis, in which the appropriate integral con-ditions describing the trim requirements were presented. The continuation analysis was first applied to a sim-ple rotor system incorporating rigid blade with flap and lag degrees of freedom with root springs to model the

stiffness of hinges. The continuation results were com-pared to published blade flap-lag stability data and ex-plained in the light of bifurcation theory. One of the key advantages of the method was the ability to identify the types of bifurcations which led to blade instabilities. This in turn helped in predicting the new secondary pe-riodic branches, leading to a better understanding of the global blade dynamics. The bifurcations discussed for this model were torus and period doubling bifurcations.

The continuation and bifurcation analysis was also applied to investigate the aeroelastic stability of a flexi-ble rotor blade, where in this case modal representation of the blade flexibility is used. This representation is widely adopted in the helicopter industry. The continua-tion results were also explained in the light of bifurcacontinua-tion theory for both forward flight and hover conditions. The results illustrated that even in the hover case the global nonlinear dynamics of the blade are quite complex. The bifurcations discussed in this case were Hopf, fold and torus bifurcations.

The analysis presented in this paper not only illus-trated that continuation and bifurcation methods are ap-plicable to studying the blade aeroelastic stability, but also confirmed that the dynamics of the blade behaviour are very complex and nonlinear, even in the hover con-dition. Although most of the complex behaviour of the blade was found at a non-operational parameter range, the continuation and bifurcation tools were essential in uncovering the global blade dynamics when multiple so-lutions coexist. Therefore, these tools offer considerable advantages in aeroelastic stability analyses of future ro-tor configurations, in particular, when the nonlinearity is significant. They are likely to offer valuable efficiencies for analyses of, for example, blades incorporating active trailing edge flaps or semi-active lag dampers.

A

The authors would like to thank AgustaWestland for funding and supporting of the work presented in this pa-per.

References

[1] A. Gessow and A. D. Crim. A method of studying the transient blade-flapping behaviour of lifting rotors at extreme operating conditions. NACA.3366, January 1955.

[2] G. Horvay. Rotor blade flapping motion. Quart. Appl. Math, Vol. 5, No. 2, 1947, pp. 149-167.

[3] G. Horvay and S. W. Yuan. Stability of rotor blade flapping motion when the hinges are tilted: Generalization of the ‘rect-angular ripple’ method of solution. Journal of the Aeronautical

Sciences, 1947, pp. 583-593.

[4] A. G. Shutler and J. P. Jones. The stability of rotor blade flap-ping motion. Aeronautical Research Council, R & M No. 3178, May 1958.

[5] O. J. Lowis. The stability of rotor blade flapping motion at high tip speed ratios. Aeronautical Research Council, R & M No. 3544, January 1963.

[6] F-S. Wei and D.A. Peters. Lag damping in autorotation by a pertubation method. Proceedings of the American Helicopter

Society 34th_{Annual Forum, pages 78–25, May 15-17, 1978.} [7] D. A. Peters. Flap-lag stability of helicopter rotor blades in

for-ward flight. Journal of American Helicopter Society, Vol. 20, No. 4, October 1975.

[8] D. A. Peters and K. Hohenemser. Application of the floquet transition matrix to problems of lifting rotor stability. Journal

of American Helicopter Society, Vol. 16, No. 2, October 1971.

[9] D.A. Peters and K.H. Hohenemser. Application of the Floquet transition matrix to problems of lifting rotor stability. Journal

of the American Helicopter Society, 26:25–33, 1970.

[10] P. P. Friedmann and D. H. Hodges. Rotary wing aeroelastic-ity - a historical perspective. Journal of Aircraft, Vol. 40, No. 6:1019–1046, November-December 2003.

[11] P. P. Friedmann and D. H. Hodges. Rotary-wing aeroelasticity: Current status and future trends. AIAA Aircraft, Vol. 42, No. 10:1953–1972, October 2004.

[12] R. L. Bielawa. Rotary Wing Structural Dynamics and

Aeroelas-ticity. AIAA (American Institute of Aeronautics & Ast, 2006.

[13] I. Chopra. Perspectives in aeromechanical stability of helicopter rotors. Vertica, Vol. 14, No. 4, 1990, pp. 457-508.

[14] D. Rezgui, P. C. Bunniss, and M. H. Lowenberg. The stability of rotor blade flapping motion in autorotation using bifurcation and continuation analysis. Proceedings of 32nd European

Ro-torcraft Forum, 2006.

[15] D. Rezgui, M. H. Lowenberg, and P. C. Bunniss. Experimental and numerical analysis of the stability of an autogyro teetering rotor. Proceedings of the American Helicopter Society 64th An-nual Forum, April 2008.

[16] D. Rezgui, M. H. Lowenberg, and P. C. Bunniss. A combined numerical/experimental continuation approach applied to non-linear rotor dynamics. Proceedings of the European

Consor-tium for Mathematics in Industry (ECMI), 2008.

[17] M. H. Lowenberg, D. Rezgui, and P. C. Bunniss. Experimental evaluation of numerical continuation and bifurcation methods applied to autogyro rotor blade aeromechanical stability.

Pro-ceedings of the ASME 2009 International Design Engineering Technical Conferences & Computers and Information in Engi-neering Conference, August-September 2009.

[18] Y.A. Kuznetsov. Elements of Applied Bifurcation Theory.

Springer-Verlag, 1995.

[19] J. Guckenheimer and P. Holmes. Nonlinear Oscillations,

Dy-namical Systems and Bifurcations of Vector Fields. Springer,

1993.

[20] E. J. Doedel, B. E. Oldeman, A. R. Champneys, F. Dercole, T. F. Fairgrieve, Yuri A. Kuznetsov, R. C. Paffenroth, B. Sand-stede, X. J. Wang, and C. Zhang. AUTO-07P : Continuation and bifurcation software for ordinary differential equations. Tech-nical report, Concordia University, Montreal, Canada, 2009. http://indy.cs.concordia.ca/auto/.

[21] Rotor Systems Group Report. Users specification for R150. Technical Report RSG/00/0217, Issue 2, AgustaWestland, May 2002.

[22] D.M. Pitt and D.A. Peters. Theoretical prediction of dynamic inflow derivatives. Vertica, Vol. 5, 1981, pp. 21-34.

[23] D. A. Peters and N. HaQuang. Dynamic inflow for practical applications. Journal of American Helicopter Society, Vol. 33, No. 4, 1988, pp. 64-68.

[24] G.H. Gaonkar and D.A. Peters. Review of dynamic inflow mod-elling for rotorcraft flight dynamics. Vertica, Vol. 12, No. 3, 1988, pp. 213-242.