Stability Analysis of Nonlinear Rotating Systems Using Lyapunov Characteristic Exponents Estimated From Multibody Dynamics

The use of Lyapunov Characteristic Exponents to assess the stability of nonlinear, time-dependent mechanical systems is discussed. Speciﬁc attention is dedicated to methods capable of estimating the largest exponent without requiring the Jacobian matrix of the problem, which can be applied to time histories resulting from simulations performed with existing multibody solvers. Helicopter ground resonance is analyzed as the reference application. Improvements over the available literature are: the problem is formulated in physical coordinates, without eliminating periodicity through multiblade coordinates; the rotation of the blades is not linearized; the problem is modeled considering absolute positions and orientations of parts. The dynamic instability that arises at some angular velocities when the isotropy of the rotor is broken (e.g., caused by the failure of one lead-lag damper, a design test condition) is observed to evolve into a large amplitude limit cycle, where the usual Floquet-Lyapunov analysis of the linearized time-periodic simply predicts instability.


Introduction
Stability assessment is a fundamental aspect of the analysis and design of any kind of dynamical system.The foundations of modern stability theory lie in Aleksandr M. Lyapunov's work [1].When considering Linear, Time-Invariant (LTI) problems, in the form ẋ = Ax (1) stability is no longer a local property of a specific solution, x(t), resulting from a specific set of initial conditions, x(t 0 ) = x 0 , but rather becomes a characteristic of the entire system.In many applications, for example those associated with rotorcraft dynamics, and specifically in helicopter rotor aeromechanics, but also in the dynamics analysis of wind [2,3] and gas turbines [4], and rotors in general [5], problems often need to be formulated as time-periodic, although often still linear, through linearization about a periodic orbit, resulting in Linear, Time-Periodic (LTP) problems.In this case, stability assessment may exploit the periodicity of the motion through the well-known Floquet-Lyapunov method [6] (see for example [7]).Asymptotic stability is associated with contraction, over a period T , of the so-called "monodromy matrix", which corresponds to the State Transition Matrix (STM) of the problem over one period.This approach is relatively popular in the rotorcraft aeromechanics community, thanks to the work of Parkus [8], who first used Floquet-Lyapunov analysis in blade flapping aeromechanics.However, it was only thanks to the advent of digital computers that the method could be applied in practice.The first attempts were proposed by Peters and Hohenemser [9], followed by Biggers [10] and Friedmann and Silverthorn [11] concerning the stability of flapping dynamics, where periodicity arises from non-axial flow, and by Hammond [12] concerning ground resonance, where periodicity may result from loss of isotropy due to mechanical failure, to name a few.A recent development by some of the authors of the present work addressed the sensitivity study of the characteristic solutions [13].
In general, however, when problems are non-linear and subjected to non-(strictly) periodic time dependence, as may occur in many transient-related problems, the ability to

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e d
evaluate the stability of reference trajectories can be extremely useful.The use of Lyapunov Characteristic Exponents (LCE) or, in short, Lyapunov Exponents (LE), has been recently proposed in the field of rotorcraft aeromechanics [14][15][16], also with a focus on their sensitivity study [17][18][19][20], and other aerospace-related applications [18].In this work, the possibility to estimate the largest LCE from time series is leveraged to support its estimation from general-purpose multibody analysis, usually formulated as a set of differential-algebraic equations (DAE), an increasingly popular tool to investigate the dynamics of complex mechanisms and systems.This would make the stability analysis of complex problems, typical for example of rotorcraft aeromechanics, feasible at low computational cost and with minimal implementation effort, if any.
The essential theoretical aspects of LCE estimation are presented in Section 2.
Section 3 briefly describes a Jacobian-based method, to support the discussion of critical aspects of LCE estimation, and the Jacobian-less method that will be subsequently used.
Section 4 illustrates the reference problem of helicopter ground resonance, using the simplified model proposed by Hammond in 1974 [12].It is selected because it represents a benchmark for the study of isotropic and non-isotropic rotors and rotor supports.Its equations of motion are formulated without any linearization.A multibody model is also formulated using a general-purpose solver by combining rigid bodies, kinematic constraints, and deformable structural components for linear springs and dampers.
Stability results from simplified and linearized models of this problem are compared with corresponding results obtained from the proposed models, and discussed in Section 5.

Theoretical Foundations
As discussed for example in [16] and references therein (e.g., [21]), LCEs indicate the rate of expansion or contraction of perturbations of a generic solution of the nonlinear differential problem with x ∈ R n , f : R n+1 → R n , and t ∈ R, along independent directions in the state space.As such, they describe the stability of the reference solution with respect to such directions.Consider a solution x(t) of Eq. (2) for t ≥ t 0 (some authors refer to it as the 'fiducial trajectory'), and a solution i x(t) of the problem where f /x x(t),t is the partial derivative of f with respect to x, evaluated along the trajectory x(t) at time t, for a perturbation i x 0 of arbitrary magnitude and direction.LCEs are defined as Each λ i is calculated from one of n linearly independent i x 0 , the equivalent of the principal directions of a LTI problem.Since Eq. ( 3) is linear time-dependent, its solution can be expressed through the State Transition Matrix (STM) so the contractive or expansive characteristic of i x is actually a property of the STM.
The STM Y(t,t 0 ) is the solution of Ẏ = f /x x(t),t Y, with Y(t 0 ,t 0 ) = I, namely the matrix that represents the evolution of the problem's state from time t 0 to t.Then Eq. ( 4) is equivalent to When all LCEs are negative, the solution is exponentially stable.When at least one LCE is positive, the solution is unstable, or leads to a chaotic attractor.When the largest LCE is zero, or the largest LCEs are zero, a limit cycle oscillation (LCO) is expected; i.e., there exists one direction, or multiple independent directions, in the state space along which the solution neither expands nor contracts.In case of multiple largest LCEs equal to zero, a higher order periodic or quasi-periodic attractor exists, e.g. a torus.
Note the analogy with the LTI case, since Y(t,t 0 ) LTI ≡ e A(t−t 0 ) , with A ∈ R n×n , and thus λ i LTI = Re(eig(A)), and with the LTP one, in which with P(t) ∈ R n×n and P(t + T ) ≡ P(t) ∀t ∈ R, where T ∈ R is the constant period, and B ∈ R n×n constant, yielding with n ∈ N. In this sense, one may consider the LCEs as sort of the eigenvalues of matrix f /x , appropriately averaged over time.For this reason, LCEs are often called the 'spectrum' of the associated problem [21], much like the spectrum of LTI problems is represented by the eigenvalues of matrix A = f /x .In most practical cases, the definition of Eq. ( 6) cannot be used in practice, because usually for t → +∞ at least

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e d
some of the elements of the STM contract to zero, in case of asymptotic stability of the solution, or expand to infinity, in case of instability, resulting in either under-or overflowing, or both.As proposed in the work of Benettin et al. [22], to estimate LCEs in practice one needs to exploit the re-orthogonalization of the local directions of evolution of the solution.Alternative approaches based on well known orthogonal decompositions (the Singular Value Decomposition (SVD) and the QR factorization, respectively [23]) have been proposed: where svd(M) indicates the singular values of matrix M, which are non-negative by definition, whereas qr(M) indicates the diagonal elements of the upper-triangular matrix R that results from the QR factorization of M (an algorithm that factors a generic real matrix M ∈ R n×n into the product of an orthonormal matrix, Q ∈ R n×n , and an upper triangular matrix, R ∈ R n×n , namely M = QR), using its realization that guarantees they are non-negative.Efficient algorithms are available for LCE estimation; for example, the continuous SVD and QR methods, and the discrete QR method [21,22,24].It is worth recalling that these methods do not easily deal with LCEs of multiplicity greater than one.This case, unfortunately frequent in mechanics, corresponds to subcritically damped oscillatory systems in LTI problems, whose complex conjugated eigenvalues share the same real part.As discussed in Ref. 25, the computability of an LCE requires that the problem possesses the property of exponential dichotomy, which in turn requires an appropriate separation between LCEs.This issue has been discussed in [26], where the real Schur decomposition [23] is proposed to effectively detect LCEs with multiplicity greater than one, although no algorithm capable of dealing with over-or underflow of the STM has been devised yet.
The possibility to extend the approach to systems of DAE, as outlined for example in [27][28][29][30][31][32], represents a promising development, in view of their use in the formulation of modern multibody dynamics based on the redundant coordinate approach.
The use of Jacobian-less methods applied to time histories resulting from multibody dynamics is very promising, as proposed by the authors in [33,34], since it intrinsically overcomes the DAE nature of the problem's formulation, and there is no need to access the Jacobian matrix of the problem, nor to evaluate it numerically, making the approach easily applicable to existing solvers.

Practical LCE Estimation
This section briefly presents the methods used in this work to estimate the LCEs or, in some cases, only the largest LCE, from problems of interest analyzed with the multibody formalism.

Discrete QR Method
The discrete QR is one of the most popular methods for the estimation of LCEs.It is based on incrementally updating the LCE estimates with the diagonal elements of the uppertriangular matrix R obtained from the QR factorization of the STM between two consecutive time steps.
Given the STM Y(t,t j−1 ) from time t j−1 to an arbitrary time t, set Y j = Y(t j ,t j−1 ).Consider then the QR factorization [23] of This way, Y j Q j−1 R Π j−1 can be used to construct the QR factorization of the STM from t 0 to t j as Y(t j ,t 0 ) = Q j R Π j by only considering incremental QR factorizations over with limited contraction/expansion in each matrix R k .Recalling that the product of two upper triangular matrices A and B, namely C = AB, is an upper triangular matrix as well, with diagonal elements c ii = a ii b ii , the LCEs can be estimated from R Π j as where j ∈ N and r ii (t j ) are the diagonal elements of matrix R(t j ) = R Π j .Thus the logarithm of each element can be incrementally computed as which helps preventing overflow or underflow in numerical computations, leading to Care must be taken to ensure that the diagonal elements of R j are positive, which can be satisfied by changing the sign of each row of the R matrix with negative diagonal element, and the corresponding column of the Q matrix; in pseudocode:

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e d
The evaluation of the STM requires matrix f /x .Alternative approaches, based on its numerical evaluation through finite differences are available, embedding it into the evaluation of the STM across each time step.In this work, the method proposed by Dieci [35] has been used in application to what has been termed the ad-hoc implementation of the problem, namely its direct formulation as a system of Ordinary Differential Equations (ODE).Matrix f /x was numerically evaluated using the complex algebra approach proposed by Squire and Trapp [36].

Jacobian-Less Methods: Maximum LCE from Time Series
The largest LCE, or Maximum LCE (MLCE) can be estimated from a time series.Among the algorithms proposed in the literature (see for example [37]), the one proposed by Rosenstein et al. [38] is used in this work.
The trajectory matrix, X, is constructed from an N-point time series x i , i = 1, ..., N, using the time delay method.Each row of matrix X is a phase-space vector, namely with with k = 1, ..., m.Thus, X ∈ R M×m , with m, M, J, and N related by where m is the embedding dimension, estimated through the False Nearest Neighbor (FNN) algorithm [39], N the length of the time series, and J represents the so-called reconstruction delay, obtained by estimating the Average Mutual Information (AMI) [40].
After constructing the trajectory matrix, the algorithm locates the nearest neighbor, X ĵ, of each point on the trajectory.It is found by searching for the point that minimizes the distance from each particular reference point, X j .The distance is expressed as where d j (0) is the initial distance from the jth point to its nearest neighbor, and • denotes the Euclidean norm.
An additional constraint is that nearest neighbors have a temporal separation greater than the mean period, T , of the time series, namely the reciprocal of the mean frequency of the power spectrum, although it can be expected that any comparable estimate, e.g., using the median frequency of the magnitude spectrum, yields equivalent results.In all the proposed cases, it is estimated as using Matlab's 'meanfreq' function, where f s is the sampling frequency of the time series.The nearest neighbors constraint is Thanks to this, each pair of neighbors can be considered as nearby initial conditions for different trajectories.
The largest Lyapunov exponent is then estimated as the mean rate of separation of the nearest neighbors.The jth pair of nearest neighbors diverge approximately at a rate given by the largest Lyapunov exponent: where C j is the initial separation.By taking the logarithm of both sides which represents a set of approximately parallel lines, for j = 1, ..., M, each with a slope roughly proportional to λ 1 .
The largest Lyapunov exponent is calculated using a leastsquares fit to the "average" line defined by where • denotes the average over all values of j.
In this work, the Jacobian-less method proposed by Rosenstein et al. is applied to time histories resulting from the simulation of the ad-hoc formulation of the problem, as an alternative to the discrete QR method, and to time histories resulting from the use of MBDyn 1 , a free generalpurpose multibody solver [41].In the latter case, this method is of particular interest, since it does not require the capability of the solver to expose the Jacobian matrix of the problem, and can thus be applied to any solver.

Example Problem: Ground Resonance
The problem of helicopter ground resonance is considered to exemplify the estimation of LCEs from a non-trivial dynamical problem, since it is characterized by non-linear, time-varying dynamics, with nonlinearities both intrinsic in the geometry of the blades' motion and possibly in the constitutive properties of the lead-lag dampers.
The analysis of the dynamic stability of a linearized model of an isotropic rotor is standard, since it results in a LTI problem when so-called Multi-Blade Coordinates (MBC) [7] are considered.When isotropy or symmetry is lost, for example when a lead-lag damper fails, the problem becomes time-periodic.In this case, its stability can be analyzed using the Floquet-Lyapunov approach.In classical applications of rotorcraft dynamics, the problem is linearized, assuming small blade lag angles, and analyzed using Floquet's theory.
In this work, however, to better highlight the ability of the proposed methods to estimate the LCEs, the problem is formulated without any small angle approximation for the rotation of the blades, and thus without a priori linearization.Furthermore, no MBC are used, an easier and more straightforward way to model a complex system using the multibody approach.

Equations of Motion
A classical model, which became a de-facto benchmark, was proposed by Hammond in his 1974 paper [12].A sketch is presented in Fig. 1.The equations of motion of the system, formulated ad-hoc, are thus for each blade, where f i ( ζi ) is the blade damping moment, with f i ( ζi ) = −c i ζi when the linear damper of [12] is considered, and with ρ = S b /m b , for the airframe.The numerical data proposed in [12] are reported in Table 1.For the airframe, they refer to quantities reduced to the equivalent motion of the hub, x h and y h , which explains the different mass values for the two coordinates.Table 1.Hammond model's data [12].

Multibody Model
The problem has been also modeled using MBDyn.The model is exactly equivalent to the simplified ground resonance problem originally proposed by Hammond.However, owing to the peculiar characteristic of the problem of having different airframe equivalent inertia terms associated with motion in the x and y directions, its modeling in MBDyn has been obtained by splitting the airframe body in two parts: • the first part is connected to the ground by a constraint that only allows its absolute displacement in the x direction; • the second part is connected to the first one by a constraint that only allows its displacement relative to the first one in the y direction; • the mass of the second part is m y , whereas that of the

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e d
Fig. 2. Sketch of the MBDyn model of Hammond's system [12].
first one is m x − m y , such that the overall mass associated with the absolute motion of the hub center in the x direction, involving the combined motion of both parts, is m x ; • the first part is connected to the ground by a spring and a damper, of characteristics k x and c x ; • another spring and damper, of characteristics k y and c y , connect the second to the first part; The rotor hub is modeled as a third, massless part, whose relative motion with respect to the second part of the airframe is a prescribed rotation about axis z with constant angular velocity; The blades are described as rigid bodies through their absolute displacement and orientation, constrained to the hub by revolute joints that only allow their relative rotation about the lead-lag hinge, whose axis is parallel to the global z axis, and thus to the axis of rotation of the rotor.Such rotation is restrained by an angular damper that represents the equivalent lead-lag damper torque.
A sketch of the model is shown in Fig. 2. The MBDyn model of Hammond's problem is available from the project's website 2 .

Isotropic System
The rotor is isotropic when all blades are identical in terms of geometry and constitutive properties.In the case at hand, the geometry and the inertia of the blades are assumed to be identical.The rotor is thus isotropic when also all the damper characteristics c i are identical, for i = 1 . . .N b .
This section compares the stability results obtained with the following models: a) the conventional, linearized time-invariant model obtained using the MBC, which acts as a reference ('Eige-  The values indicated as "MLCE MBDyn" refer to the maximum LCE estimated from the time series computed by the multibody model, whereas those indicated as "LCE QR ad-hoc" have been obtained from the ad-hoc formulation using the QR algorithm.All results are negative, which indicates asymptotic stability.These results are in line with those proposed by Hammond in [12], which were obtained from a linearized model using MBC. It is worth noticing the set of values at about −1.875 s −1 from 65 rpm onwards, which actually represents two coincident eigenvalues, associated with the 'collective' (all blades rotating by the same amount) and so-called 'reactionless' (blades alternately rotating by opposite amounts) rotor MBC.Those two real parts of the corresponding eigenvalues are not affected by the angular velocity of the rotor, since the corresponding motions are mechanically decoupled from the motion of the airframe.As such, usually they are not considered when the problem is formulated using MBC; they appear in the solution presented here because it is formulated using the four physical blade rotation angles as coordinates.

Loss of Isotropy: One Damper Inoperative
Current Certification Standards (CS 27 [42] and CS 29 [43], respectively for small and large rotorcraft) require that "the rotorcraft may have no dangerous tendency to oscillate on the ground with the rotor turning" (CS 27.241/29.241),and that "the reliability of the means for preventing ground resonance must be shown either by analysis and tests, or reliable service experience, or by showing through analysis or tests that malfunction or failure of a single means will not cause ground resonance" (CS 27.663(a)/29.663(a)).
For this purpose, the case of one lead-lag damper inoperative needs to be evaluated, to check the stability margins of the rotorcraft within the desired envelope of rpm and other structural parameters (weight and center of mass position, landing gear equivalent stiffness and damping, and so on).In the present analysis, we focus on checking the suitability of the proposed methods in analyzing a problem that no longer can be made time-independent, and that according to the linearized analysis proposed in [12] based on the Floquet-Lyapunov theory results in extensive regions of instability.
In this section we compare the stability results obtained with the following models: a) the conventional, linearized time-periodic model obtained using the MBC, although their use no longer makes the problem LTI because of the loss of isotropy; this solution is used as a reference ('Floquet'); b) the results obtained using the Jacobian-based discrete QR method to estimate all the LCEs of the ad-hoc model ('LCE QR ad-hoc'); c) the results obtained using the Jacobian-less method for the estimation of the MLCE from the time series computed by the multibody and ad-hoc models ('MLCE MBDyn' and 'MLCE ad-hoc').
The latter two cases are identical to those considered for the isotropic rotor problem.
The lag damper of blade 3 is made inoperative, i.e., c 3 = 0. Figure 4 compares the real part of Floquet's exponents resulting from the linearized, time-periodic analysis, with the LCEs estimated from the nonlinear ad-hoc problem using the discrete QR method.One may observe that where Floquet's method can only find positive real parts, indicating instability of the linearized model, the largest of the LCEs estimated from the geometrically exact nonlinear problem has (almost) zero value, which is indicative of an LCO, as discussed later when the phase space trajectories of the blades are presented in Fig. 8 for an unstable case.Indeed, the blade can only rotate by a finite, although large amount about the lag hinge; thus, a motion limited in amplitude results instead of the indefinitely divergent motion predicted by the analysis   of the linearized periodic model.
Figure 5 compares the largest real part of the Floquet exponent with the Maximum LCE, estimated from the time histories obtained by integrating the multibody problem using the Jacobian-less method.Identical results have been obtained with the ad-hoc formulation of the problem, as one would expect, considering that the resulting time histories are quite similar; they only slightly differ because different numerical integration schemes have been used.For the ad-hoc formulation, Matlab's ode23() implementation of the explicit Runge-Kutta scheme proposed by Bogacki and Shampine [44] was used.With MBDyn, a second-order accurate, implicit A-stable multistep scheme with tunable algorithmic dissipation and asymptotic spectral radius ρ ∞ = 0.6 [45] was used.
Figure 6 is analogous to Fig. 5; the curves differ in that the analyses used to estimate the MLCE were executed for a fixed time duration.Consequently, when the MLCE is negative and too large in modulus, the analysis results quickly underflow, making the estimation less accurate.For this reason, particular care is required in defining the appropriate duration of the analysis.
Figures 7 and 8 show the phase space of the 6 coordinates of the ad-hoc problem at two specific values of rpm and initial configuration.Figure 7 refers to an asymptotically stable case, with Ω = 124.2rpm and an initial solution of x 0 = 0.1 m for the airframe displacement in the x direction and everything else zero.The plot of blade 3, the one missing the damper, shows larger amplitude and much slower decay of the oscillations compared with the others.It is worth noticing that the phase space portraits of blades 2 and 4 are  identical if the signs of the axes are exchanged.The circles indicate samplings spaced by one period, i.e., by one rotor revolution.
Figure 8 refers to an unstable case, with Ω = 242.4rpm and the same initial condition.As illustrated by the (nearly) zero MLCE of Fig. 5, the solution evolves towards a LCO, highlighted by the circles obtained by sampling the solution at each revolution that describe a periodic orbit.The red and blue spots in the x and y plots correspond to the solution entirely spanning that portion of the phase space, eventually resulting in the aforementioned periodic orbit.The ampli-tude of the motion of blade 3, the one with the inoperative damper, is of the order of ±π/2 (−1.9 radian to 1.5 radian); quite large and surely destructive for a helicopter in terms of induced structural loads, yet undoubtedly limited.
Figure 9 shows the amplitude of blade 3 rotation with the failed damper as a function of the rotor angular velocity.The unstable rpm region can be well appreciated.As the rpm increases from zero, when unstable conditions are reached at about 210 rpm the oscillation immediately jumps to a very large extension.Then, as the rpm increases further, it slowly reduces until stable conditions are reached again  slightly above 300 rpm.It is worth stressing that the abnormal amplitude motion experienced in the region of rpm that formally corresponds to a LCO identifies an operation range that cannot be physically acceptable, as the resulting dynamic loads would likely cause the breakdown of any structure.
In conclusion, from a design point of view the results obtainable with the linearized model and Floquet-Lyapunov's analysis suffice in predicting instability.However, the correct identification of a LCO-type solution provided by the proposed analysis gives further insight into the characteristics of the phenomenon.

Nonlinear Constitutive Law
The case of a more realistic damper model, characterized by a nonlinear constitutive law, is considered.The relatively simple model used for example in [16,46], originally proposed by [47,48], is here considered, namely , with a moment-saturation limit angular velocity ζL = 1 deg•s −1 , and C L = c i .Namely, the damper moment, shown in Fig. 10, is second-order polynomial for angular rates ζi that do not exceed, in norm, a threshold ζL , and otherwise constant.The parametric stability of the isotropic problem at fixed rpm is then studied by varying the linear contribution in Eq. ( 26), namely the term −C L ζi , with C L ∈ [0, c i ], as proposed in [16], after the insurgence of a LCO for C L = 0 was observed in [46].The results presented in [16] in terms of LCEs are here reproduced and presented in Fig. 11, considering the ad-hoc model, which includes the geometric nonlinearities associated with the finite motion of the parts, and the multibody model, which on top of that adds the formulation of the blade motion with respect to the absolute reference frame, in of a set of DAEs.The two analyses show essentially identical results, since the resulting time histories are quite similar.The stability of the trivial solution, namely zero blade leadlag rotation and zero airframe displacement, is evaluated.The largest LCE is zero for values of the linear contribution to the characteristic moment of the damper, −C L ζ, that range from zero to about 35% of the nominal value, indicating the existence of a limit cycle oscillation that is confirmed also in case of complete geometric nonlinearity in the kinematics of the blades.From that point on, the LCE becomes progressively negative with steady slope, reaching a value of about -1 s −1 when the linear contribution −C L ζ reaches its nominal value.
Figure 12 shows the sensitivity of the LCE estimation to the amplitude of the initial angular velocity perturbation of blade 3, ζ3 (0), for C L = c i (the nominal value).Below a magnitude of about 0.6 rad/s, the MLCE is negative for all values of rpm.For higher initial values, a region appears, the red, flat area that corresponds to an MLCE value of about 0, where the aforementioned LCO occurs.The region grows, extending over a broader rpm range, for initial values of ζ3 (0) of increasingly larger magnitude.

Conclusions
This paper discussed the use of estimated Lyapunov Characteristic Exponents to evaluate the stability of solutions of complex mechanical systems.Its elements of novelty are: • the Maximum LCE can be estimated from generalpurpose multibody formulations, without the need to modify existing solvers, by resorting to Jacobian-less methods that make use of time series only; • very effective algorithms based on the discrete QR factorization and capable of estimating all LCEs can be im- plemented without the explicit Jacobian matrix in analytical form, by resorting to effective and accurate methods based on finite differences; • the problem of helicopter ground-resonance, formulated without a priori linearization of the blade motion or equivalently by general-purpose multibody dynamics, results in a limit-cycle oscillation, whereas the conventional linear, time-periodic analysis simply indicates instability; • the presence of a limit-cycle instability when considering a certain type of nonlinear lead-lag damper presented in the literature is confirmed also when the motion of the blades is not linearized.

Figure 3
Figure 3 illustrates how the stability indicators of the isotropic rotor change with the angular velocity of the rotor.The values indicated as 'Eigenanalysis' refer to the real part of the eigenvalues of the monodromy matrix obtained from a linearization of the problem about the trivial solution of zero blade rotation.The corresponding results obtained using the discrete QR algorithm from the ad-hoc problem yield identical results.Indeed, being the reference solution of zero blade motion always asymptotically stable, the absence of a priori linearization is inessential.The values indicated as "MLCE MBDyn" refer to the maximum LCE estimated from the time series computed by the multibody model, whereas those indicated as "LCE QR ad-hoc" have been obtained from the ad-hoc formulation using the QR algorithm.All results are negative, which indicates asymptotic stability.These results are in line with those proposed by Hammond in[12], which were obtained from a linearized model using MBC.It is worth noticing the set of values at about −1.875 s −1 from 65 rpm onwards, which actually represents two coincident eigenvalues, associated with the 'collective' (all blades rotating by the same amount) and so-called 'reactionless' (blades alternately rotating by opposite amounts) rotor MBC.Those two real parts of the corresponding eigenvalues are not affected by the angular velocity of the rotor, since the corresponding motions are mechanically decoupled from the motion of the airframe.As such, usually they are not considered when the problem is formulated using MBC; they appear in the solution presented here because it is formulated using the four physical blade rotation angles as coordinates.

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e dFig. 7 .
Fig. 7. Phase space of ad-hoc solution with one damper inoperative, in asymptotically stable conditions (Ω = 125 rpm).

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e dFig. 9 .
Fig. 9. Amplitude of blade 3 limit cycle oscillation with failed damper.

A c c e p t e d M a n u s c r i p t N o t C o p y e d i t e dFig. 12 .
Fig. 12. MLCE of isotropic case with nonlinear ad-hoc model and nonlinear damper, sensitivity to amplitude of initial perturbation.