An amplitude phase formulation for nonlinear modes and limit cycles analysis through invariant manifolds

The aim of this paper is to show how the concept of nonlinear modes can be used to characterize periodic orbits and limit cycles in multi-degree-of-freedom nonlinear mechanical systems. In line with previous studies by Shaw and Pierre, the concept of nonlinear modes is introduced here in the framework of invariant manifold theory for dynamical systems. A nonlinear mode is deﬁned in terms of amplitude, phase, frequency, damping coefﬁcient and mode shape, where the last three quantities are amplitude and phase dependent. An amplitude-phase transformation is performed on the nonlinear dynamical system, giving the time evolution of the nonlinear mode motion via the two ﬁrst-order differential equations governing the amplitude and phase variables, as well as the geometry of the invariant manifold. The system of formulation adopted here is suitable for use with a Galerkin-based computational procedure. The existence and stability of periodic orbits such as limit cycles on the associated invariant manifolds can be studied from the differential equations governing the amplitude and phase variables. The examples given here involve adding gyroscopic and/or ‘‘negative’’ nonlinear damping terms of Van der Pol type, and nonlinear restoring force to the system equations.


Introduction
Nonlinear modes (NMs) are efficient tools for analyzing the behavior of vibrating mechanical systems [1][2][3][4][5][6].In a series of papers by Shaw and Pierre [5][6][7][8], the NM concept was related to invariant manifold theory for dynamical systems.A NM in a (regular enough) n-dimensional, autonomous, second-order nonlinear mechanical system is any motion which evolves in a two-dimensional invariant set.This set contains the equilibrium point (zero oscillation) and is tangent to the plane-eigenspace corresponding to the linear mode of the associated linearized system (small oscillations).The invariant manifold geometry is characterized by a set of partial differential equations (PDEs) with respect to the ''master coordinates'' (all displacements and velocities are related to a single pair of displacement and velocity terms, called the ''master coordinates'').The time evolution of a NM motion is given by the remaining second-order nonlinear oscillator governing the master coordinates.The amplitude-phase transformation performed on the master equation, based on the linear modal frequency in Ref. [8], takes computational advantage of the 2p-periodic character of the parametric expressions for the invariant manifold in terms of the phase variable.
We recently proposed [9] an amplitude-phase formulation for characterizing a NM in conservative nonlinear systems in the line with the approach developed in Ref. [3].As in the linear case, a NM is described in terms of the amplitude, mode shape and frequency, where in this study the distinctive feature is the fact that the last two quantities are amplitude and phase dependent.The time evolution of the NM motion was then given by a first-order differential equation governing the total phase motion, from which the period of the oscillation can easily be deduced.It was established that the frequency and mode shape functions solve a 2p-periodic (with respect to the phase variable) nonlinear differential eigenvalue problem.This method also gives a parametric description of the associated invariant manifold.
In the present study this formulation is extended to autonomous mechanical systems including displacement and velocity nonlinear terms.The amplitude-phase transformation involves a new function to characterize the amplitude modulation, that we will call a ''damping function'' (in the wide sense).This function also depends on the amplitude and phase variables.The dynamics is now described by two first-order nonlinear differential equations, the one governing the total phase motion, and the other the amplitude motion.It is observed that in some cases periodic solutions may exist for the amplitude equation, denoting the presence of stable or unstable ''limit cycles'' on the invariant manifold.The procedure is first discussed in the case of single-degree-offreedom (sdof) oscillators in Section 2. Three cases are considered: a sdof oscillator without damping but with inertial nonlinearity and a nonlinear restoring force, a sdof oscillator with a linear damping term, unit mass and a nonlinear restoring force, a sdof oscillator with ''negative nonlinear damping'' of the van der Pol type yielding two stable and one unstable limit cycles.In this last example, it is observed that a small parameter approach leads to predictive analytical results as to the behavior of the oscillator.In Section 3, the procedure for n-dimensional mechanical systems is presented.For computing the NMs, the use of a Galerkin procedure based on truncated Fourier series with respect to the phase variable, combined with polynomial functions for the amplitude variable, is also discussed.A two-dimensional example is first presented.Four configurations are considered, depending on the properties of the damping terms: a conservative non-gyroscopic system, a conservative gyroscopic system, a dissipative system with proportional damping, and a dissipative system with non-proportional damping.Lastly, a system of two coupled van der Pol equations is discussed, including the stability results of the corresponding limit cycles on the associated invariant manifolds.

Undamped case with inertial nonlinearity
Let us consider an undamped sdof oscillator with inertial nonlinearity and a nonlinear restoring force.Let T ¼ 1 2 mðqÞ _ q 2 and U ¼ UðqÞ be the kinetic energy function and the potential energy function, respectively, where q denotes the displacement.Using the notation g q ðqÞ¼ðqg=qqÞðqÞ, the Euler-Lagrange equation reads giving the equation of motion where f ðqÞ¼U q ðqÞ.It will be assumed in what follows that the functions m and f are regular enough and mðqÞXm 0 40, f ð0Þ¼0 and qf ðqÞX0.
We will now establish that a periodic orbit of Eq. ( 2) can be expressed, either Numbers a and j are two constants which set the initial conditions of the motion.Eqs. ( 3) and ( 5) introduce two amplitude-phase transformations of the equation of motion (2).In both cases the frequency function O and the bias function b need to be found.These expressions differ from the classical polar transformation, since they involve a non-constant frequency Oða; fÞ.
In expressions (3), b is a scalar function which depends only on a, and O is a positive scalar function which depends on a and f.As we shall see, O is an even 2p-periodic continuous function with respect to f and hence with Fourier series In expressions (5), the bias function b depending on a and f will be sought for, arbitrarily, along with the frequency function O, in the form of an even p-periodic function with respect to f with Fourier series In both cases, it can be seen from Eqs. ( 4) and ( 6) that the period of the motion is given by After some algebraic and differential manipulations, substituting Eq. (3) into Eq.( 2) gives where with the constraint equation i.e.
Eq. ( 13) defines the bias function bðaÞ and, together with the assumption qf ðqÞX0, means that O 2 , given by Eq. (11), is a positive continuous 2p-periodic function with respect to f for fixed a with Fourier series (7).The above result is well known in the case of mðqÞ¼1, see for example Ref. [13].Unfortunately, it does not seem to be applicable to multi-dof (mdof) systems.Substituting Eq. ( 5) into Eq.(2) gives d df ðmðaxðfÞÞO 2 sin 2 fÞÀ 2 a sin ff ðaxðfÞÞ where now xðfÞ¼cos f þ bða; fÞ.Eq. ( 14) looks rather like Eq. (10) but it contains an additional term which depends on the first two derivatives of the bias function b.Eq. ( 14) cannot be solved analytically and apparently we have to solve, for fixed a, one equation for two unknown functions O 2 and b.Nevertheless we look for even p-periodic functions (see Eq. ( 8)).Consequently, using the decomposition ¼ m e ðfÞþm o ðfÞ, where ðÁÞ e (respectively ðÁÞ o ) denotes the even (respectively, odd) cosine terms in Fourier series, Eq. ( 14) can be re-written as where Eq. (15) (respectively, Eq. ( 16)) stands for the odd (respectively, even) sine terms of Eq. ( 14).If m o 0 and f e 0, then bða; fÞ¼0 is a solution of Eq. ( 16), and Eq.(15) reduces to Eq. (10), yielding the same expression for Eqs.(3) and (5) and O is p-periodic with Fourier series (8).
Except for some special cases, Eqs.(15) and ( 16) have to be solved numerically and this can be done using a Galerkin procedure based on the truncated Fourier series of Eq. ( 8) as we shall see later on.The formulation ( 5) and ( 6) can be extended to the mdof case [9] and will be used in this paper.

Odd restoring force with linear damping
Consider a sdof oscillator with linear damping, unit mass and a nonlinear restoring force We will assume that the function f is regular enough and, for the sake of simplicity, that f ðqÞ¼Àf ðÀqÞ.
If ca0, the amplitude of the motion is time dependent.To take this fact into account, a solution of Eq. ( 17) is sought for in the following: Substituting Eqs. ( 18) and (19) into Eq.( 17) and equating the sine and cosine terms separately leads to the following PDEs in terms of the two variables v and f for the two unknown functions O and x If f ðqÞ¼kq, then O 2 k Àðc 2 =4Þ and x Àc=2 as was to be expected.Otherwise, it can be established that in the case of positive (respectively, negative) damping, c40 (respectively, co0), we have xp0 (respectively, xX0Þ so that the zero solution is asymptotically stable (respectively, unstable).
Except for some special cases, Eqs. ( 21) and ( 22) have to be solved numerically, and this can be done using a Galerkin procedure based on the truncated Fourier series of Eq. (20).
The transformation Eqs. ( 18) and ( 19) can be applied to systems with nonlinear damping as well as to systems including auto-oscillations.In this last case, ''limit cycles'' can be studied as we shall see in the next section.

Application to limit cycles
From Eq. ( 19), it follows that where tðv; fÞ¼ðxðv; fÞ=Oðv; fÞÞ can be viewed as a ''generalized damping rate function''.Since the right-hand side of Eq. ( 23) is p-periodic with respect to f, periodic solutions v Ã ðfÞ¼v Ã ðf þ pÞ may exist with some x and O (one necessary condition being that tðv; fÞ does not keep a constant sign).It follows that qðtÞ¼v Ã ðfðtÞÞ cos fðtÞ will be a T-periodic function with period The stability analysis of this periodic motion (''limit cycle'') will be easily deduced from the variational equation associated with Eq. (23).
Example.Let us consider the nonlinear sdof oscillator of the van der Pol type where is a positive (small) parameter.
Substituting Eqs. ( 18) and (19) into Eq.( 25) gives similar equations to Eqs. ( 21) and ( 22) 1 2 where pðqÞ¼Àð1 À 22 3 q 2 þ 8q 4 À 32 15 q 6 Þ.One can check that solutions of Eqs. ( 26) and ( 27) are such that It follows from Eq. ( 26), that and, from Eq. ( 27), we deduce that with 28) and (29) into Eq.(23), the existence of limit cycles for Eq. ( 25) can be deduced from the existence of equilibrium points of the averaged equation ( [11], recalling that x ¼ e x and hence t ¼ et) The amplitudes of the limit cycles are given by the three positive roots of the equation Þ=qvÞ40), and the third root v Ã 3 ¼ ffiffi ffi 3 p gives a stable limit cycle (ðqx 0 ð ffiffi ffi 3 p Þ=qvÞo0).The phase picture of the stable limit cycles obtained by directly simulating Eq. ( 25) is shown in Fig. 1.From (24), the period of the limit cycles is given by Table 1 compares, for various values of the parameter , the period of the limit cycle v Ã 1 ¼ 1 with the ''exact'' period of the limit cycle obtained by directly simulating Eq. ( 25) with a shooting method.

mdof Systems
We now consider a mdof system of the form where QðtÞ is an n-vector, ½M is a non-singular symmetric n Â n-matrix and F is a (sufficiently regular) vector function with dimension n such that Fð0; 0Þ¼0.We want to characterize the NMs in the framework of the invariant manifold theory [5], using an amplitude-phase transformation as for sdof oscillators.

Undamped case
In a previous work [9] we considered the case in which the nonlinear force function F does not depend on _ Q.We briefly summarize here the main results.
Assume firstly that F is an odd function (FðQÞ¼ÀFðÀQÞ).A normal mode is defined from a family of periodic solutions to Eq. (32) sought in the form QðtÞ¼aWða; fðtÞÞ cos fðtÞ, where the scalar function f satisfies the differential equation By analogy with the linear case (where Eqs.(33) and (34) hold with constant W and O), the vector function, W and the scalar function O, will be referred to here as the modal vector and the resonance frequency of the nonlinear normal mode, respectively, and the scalar function f, will be referred to here as the total phase of the motion.The one-dimensional differential equation (34) governing the total phase motion f, will define the dynamics of the periodic response and the scalar quantities a (a40) and j (j 2½0; 2p) will set the initial conditions of the vibration in the phase space.To ensure that the parameter a appropriately characterizes the where and ðÁÞ f denotes partial differentiation with respect to f.The differential rule ðO 2 Þ f ¼ 2OO f has been used to work with the unknown function O 2 in place of O.
Eq. ( 36), which is 2p-periodic with respect to f, together with the normalization condition (35), can be viewed as a nonlinear 2p-periodic eigenvalue-eigenvector problem.Thus, the modal vector W and the frequency O will be searched for as 2p-periodic functions with respect to f for fixed a.More specifically, due to the odd symmetry of the function F, and the quadratic form (in W) of the normalization constraint (35), these functions will be searched for as even periodic functions with respect to f with period p and hence with Fourier series containing only cosine terms of even order.It follows that the Fourier series of QðtÞ in terms of fðtÞ will contain only cosines terms of odd order.It can be established that if Q T FðQÞX0 (which is weaker than the monotony assumption made in Ref. [9]) necessarily leads to O 2 X0.The period of the vibration in mode motion depends only on the amplitude a and follows easily from Eq. (34): Eq. (33) together with _ QðtÞ¼aOða; fðtÞÞðW f ða; fðtÞÞ cos fðtÞÞ À Wða; fðtÞÞ sin fðtÞÞ, (39) define a ''synchronous'' periodic oscillation in the sense of Ref. [1].The modal line in the configuration space can be either straight or curved.It should be mentioned that this formulation gives also a characterization of the NM in the framework of the invariant manifold in the phase space [5,7,8].The invariant manifold coincides with the set, in the phase space, defined by Eqs.(33) and (39) taking the initial conditions a and j to be independent variables.If FðQÞa À FðÀQÞ, due to the presence of quadratic terms for example, Eq.(33) must be replaced by [9] QðtÞ¼aðWða; fðtÞÞ cosðfðtÞÞ þ Bða; fðtÞÞÞ, where the term Bða; fÞ, such that Bða; fÞ¼Bða; ÀfÞ¼Bða; f þ pÞ; takes into account the even cosine terms in the nonlinear oscillation.
To prepare the next section, an equivalent formulation, more convenient for the numerical analysis, consists in seeking (Q; _ Q) in the form QðtÞ¼aXða; fðtÞÞ; _ QðtÞ¼aYða; fðtÞÞ; where fðtÞ is again given by Eq. (34).Substituting Eq. (41) into Eq.(32), using Eq.(34) we obtain a system of first-order nonlinear differential equations in the variables f for the unknowns X; Y; O, for which 2p-periodic solutions are planned.To complete the above equations a normalization condition must be added.We use the following result: any regular, even, 2p-periodic vector function XðfÞ can always be written XðfÞ¼X oc ðfÞþX ec ðfÞ, where X oc (resp., X ec ) stands for the odd cosine (resp., even cosine) terms in the Fourier series of XðfÞ.Moreover, there exists one and only one even and p-periodic vector function in f, W c ðfÞ, such that X oc ðfÞ¼W c ðfÞ cos f.It follows that the normalization constraint (35) can again be retained and rewritten as Finally, we seek periodic solutions of Eqs. ( 42)-(44) such that Oða; fÞ¼Oða; ÀfÞ¼Oða; f þ pÞ. (47)

The general case
Return now to the general system of Eq. ( 32).This may include systems having damped modes or systems with auto-oscillations.In both cases the amplitude is no longer constant but time-modulated.As in the sdof systems (see Section 2), we must consider a supplementary first-order equation governing the amplitude modulation.
More specifically, we focus on motion where the displacements and velocities, ðQðtÞ; _ QðtÞÞ, are related to a single pair of amplitude and phase variables, ðvðtÞ; fðtÞÞ, according to QðtÞ¼vðtÞXðvðtÞ; fðtÞÞ; _ QðtÞ¼vðtÞYðvðtÞ; fðtÞÞ; where X and Y are n-vector functions.The amplitude and phase variables are governed by the two first-order differential equations _ fðtÞ¼OðvðtÞ; fðtÞÞ; _ vðtÞ¼vðtÞxðvðtÞ; fðtÞÞ; fð0Þ¼j; vð0Þ¼a, ( where O and x are scalar functions.Numbers j (2½0; 2p) and a (40) are two given constants which set the initial conditions of the motion.As in Section 3.1, the four functions X, Y, O and x will define a modal motion for Eq.(32), provided that they are 2p-periodic in the phase variable f.Again, the two-dimensional invariant manifold associated with the NM in the phase space reads while the two first-order differential equations (49) govern the dynamics of the motion on the invariant manifold.
The two scalar functions, O (the frequency function) and x (the (generalized) damping function), play the same role as in the sdof case (see Section 2).
Substituting Eq. (48) into Eq.(32) and using Eq. ( 49) and the chain rule, yields a set of first-order nonlinear PDEs in the two variables ðv; fÞ, The PDEs (51) and (52) are independent of time and contain the two unknown vector functions X and Y, their partial derivatives with respect to v and f and the two unknown scalar functions O and x.
To characterize the four unknown functions X, Y, O and x, it is necessary to add two constraint scalar equations to Eqs. (51) and (52) (normalization conditions).Due to the 2p-periodicity in the phase variable f, the function X can always be broken up according to where ð:Þ oc (ð:Þ ec , ð:Þ os , ð:Þ es , respectively) denotes the odd cosine (even cosine, odd sine, even sine, respectively) terms in the corresponding Fourier series.Moreover there exist (one and only one) even and p-periodic vector functions W c ; W s such that X oc ðv; fÞ¼W c ðv; fÞ cos f; X os ðv; fÞ¼W s ðv; fÞ sin f.
It follows that we can adopt the normalization condition (usually set for linear systems): W T c ½MW s ¼ 0, which can be rewritten as Finally, a NM of the system (32) is obtained by solving Eqs. ( 51 which can be interpreted in terms of amplitude modulation coefficient (v), eigenvectors (W c ; W s ) and frequency (O), like for a linear system.

Computing the nonlinear modes
The challenge is now how to solve the PDEs (51)-(54) for the four unknown functions X, Y, O and x.These equations have been formally established for ðv; fÞ2R þ Â½0; 2p.In the approximation procedure, the domain over which the PDEs will be solved will be limited to ½0; v max Â½0; 2p, where v max is a positive constant.

Each one of these linear modes can be written
Xðv; fÞ¼W c cos f À W s sin f; Oðv; fÞ¼o and xðv; fÞ¼Z: where the complex frequency The eigenvector, W ¼ W c þ iW s , is next given by solving the two following equations: under the normalization constraints We focus now on calculating the NMs which coincide at small values of v (i.e.v % 0) with one of the modes of the linear system (57).
Since the unknown functions are 2p-periodic and p-periodic with respect to the variable f (see Eqs. ( 55) and ( 56)), a Galerkin discretization procedure based on the truncated Fourier series with respect to the variable f can be used advantageously here.The trigonometric functions are combined with the canonical polynomial functions with respect to the variable v. Hence, the unknown functions are expanded into X oc ðv; fÞ%X c lin ðfÞþ X os ðv; fÞ%X s lin ðfÞþ Y oc ðv; fÞ%Y c lin ðfÞþ Y os ðv; fÞ%Y s lin ðfÞþ  51) and ( 52)) and v p cos 2mf (in Eqs. ( 53) and ( 54) algebraic equations in the unknown X c p;2kþ1 , X c p;2k ; ...; x p;2k .For fixed N v and N f , a Newton-Raphson method with a zero starting point has been implemented to solve these equations.
Based on some symmetrical properties of the nonlinear function F, the number of unknowns can be significantly reduced.More specifically, if F is a polynomial vector function with only odd degree terms, the approximation procedure can be simplified, keeping in Eqs. ( 64)-(66) only terms with even degrees with respect to the variable v and setting X ec 0, X es 0, Y ec 0 and Y es 0.

Linear damping and nonlinear restoring forces
Let us consider the two-dimensional mechanical system presented in Refs.[8,9] with o 1 ¼ 0:689 and o 2 ¼ 3:244.Four configurations have been studied, depending on the properties of the damping terms: non-gyroscopic and gyroscopic conservative systems, and dissipative systems with proportional and non-proportional damping.With each configuration, the NM which coincides at small values of v with the second mode of the underlying linear system was computed by solving Eqs. ( 51)-(54) over ½0; v max Â½0; 2p with v max ¼ 1:9a s described in Section 4. Due to the symmetrical properties of the system, we have X ec 0, X es 0 imposed and only terms with even degrees with respect to the variable v have been kept in Eqs. ( 64)-(66).The order of truncation was selected, based on the normalized mean square error over ½0; T max defined by , where Q ¼ðq 1 ; q 2 Þ.Q a denotes the modal motion given by the amplitude-phase formulae (Eq.( 48)) by solving Eqs.(49) numerically with the initial condition a ¼ v max and j ¼ 0. Q e denotes the corresponding motion obtained by solving the equation of motion Eq. (67) with the initial condition Q e ð0Þ¼Q a ð0Þ and _ Q e ð0Þ¼ _ Q a ð0Þ.The range of integration T max depends on the behavior of the modal motion under investigation.

Non-gyroscopic conservative system
Parameter values: Linear mode under consideration (see Section 4): Z ¼ 0, o ¼ 3:244, w c ¼ð0; 0Þ T and w s ¼ð0; 1Þ T .The corresponding NM is shown in Figs. 2 (a)-(d).All the coefficients in the x-expansion (Eq.( 66)) as well as those in the X oc -expansion (Eq.(64)) were found to be equal to zero.An X-expansion with only cosine terms was obtained upon taking w c ¼ð0; 1Þ T and w s ¼ð0; 0Þ T .This result is in agreement with the amplitudephase approach performed in Ref. [9], where the functions x and X os were both taken to be equal to zero.Since function x is equal to zero, the modal motions are periodic and the period versus initial amplitude a (see Eq. ( 38)) is plotted in Fig. 2(d).The nonlinearity affects not only the dynamical properties of the mode (see the behavior of the period TðaÞ with respect to the amplitude level a) but also the geometrical properties of the mode (see the behavior of X with respect to the amplitude level v).The associated invariant manifold is therefore not a linear subspace (see Fig. 6(a)).
The normalized mean square errors computed with T max ¼ Tðv max Þ (given by Eq. ( 38)) with different values of N v and N f are given in Table 2. Orders N v ¼ 4 and N f ¼ 2, give sufficiently accurate approximations over ½0; 2pÂ½0; v max .
The corresponding NM is shown in Figs. 3 (a)-(d).The first component of X shows the gyroscopic effect.The modal motions are again periodic (see Fig. 6(a)).Contrary to the non-gyroscopic conservative case, the coefficients in the x-expansion (Eq.(66)) were not found to be both equal to zero, but to have small values.It was checked that imposing x 0 yields the same result.The normalized mean square error was rðTðv max ÞÞ ¼ 0:00433 with N v ¼ 4 and N f ¼ 2.
The corresponding NM is shown in Figs.4(a)-(d).As was to be expected, the damping function x was found to differ significantly from zero, but kept a constant (negative) sign starting at x ¼À0:649 at v ¼ 0 and increasing slowly with v.The corresponding invariant manifold, including the modal motion with initial Table 2 Normalized mean square error conditions vð0Þ¼v max and fð0Þ¼0, is given in Fig. 6(c).The motion decreases slowly up to the equilibrium point.The normalized mean square error was computed with T max ¼ 10Tðv max Þ, yielding rðT max Þ¼0:00962 with N v ¼ 4 and N f ¼ 2.
The corresponding NM is shown in Figs. 5 (a)-(d).The behavior of the frequency and damping functions was found to be similar to that obtained in the proportional damping case, except for the first component of the vector X.Note that with large v, the difference is smaller but still persists.The corresponding invariant manifold, including the modal motion with initial conditions vð0Þ¼v max and fð0Þ¼0, is given in Fig. 6(c).Again, the motion decreases slowly up to the equilibrium point.
The normalized mean square error was computed with T max ¼ 10Tðv max Þ yielding rðT max Þ¼0:0103 with The invariant manifolds were plotted in Figs.6(a)-(d) in the phase subspace ðQ 1 ; Q 2 ; _ Q 2 Þ using the same viewpoint and the same plot ranges.For each configuration, the invariant manifold of the linearized system is also shown on the faces of the bounding box.For conservative systems, the modal motions are periodic orbits.Each periodic orbit is characterized solely by its amplitude, and the invariant manifold can be said to be the union of all these periodic orbits.These properties can be used to characterize and compute a NM (see Ref. [12]).In the case of dissipative systems, the modal motions are orbits decreasing slowly up to the equilibrium point.Here, again the invariant manifold can be said to be the union of orbits with initial conditions vð0Þ¼v max and fð0Þ¼j 2½0; 2p.

A Van der Pol system
Let us consider the two coupled van der Pol equations Fig. 7 shows the computed frequency and damping functions.At v ¼ 0, these functions coincide with the modes of the linearized system.Both frequency functions increase with v, reflecting the hardness behavior of the model in the two NMs.Both damping functions decrease with v, beginning with positive values and taking negative values at high values of v.The solutions of Eq. ( 68) on the invariant manifolds can therefore be expected to be periodic.The dynamic behavior of the modal motions are governed by the two first-order differential equations (49).A periodic motion may occur on the invariant manifold if there exits a periodic solution to Eq. ( 23).The existence of a periodic solution is deduced from the existence of an equilibrium point in the averaged equation (using the average principle in the context of perturbation theory, see Ref. [ Mean damping rates hti versus v are plotted in Fig. 8 for the two NMs.Eq. ( 69) takes the equilibrium point v Ã 1 ¼ 1:59 in the first NM and v Ã 2 ¼ 1:83 in the second NM.Each equilibrium point characterizes an assymptotically stable (d=dvhtiðv Ã i Þo0) limit cycle on the associated invariant manifold.The approximations can be improved by solving Eq. ( 23).The harmonic balance method yields the following truncated expansion: v Ã 1 ðfÞ¼1:5957 À 0:0059 cos 2f À 0:0010 cos 4f À 0:0430 sin 2f þ 0:0004 sin 4f, v Ã 2 ðfÞ¼1:8657 þ 0:0461 cos 2f À 0:01761 cos 4f þ 0:2146 sin f À 0:0043 sin 4f.giving the limit cycle approximations (i ¼ 1; 2) The last question which now arises is whether or not the limit cycles are stable in the phase space.To answer this question Floquet's theory (see Ref. [11]) is applied.Rewriting Eq. ( 68) in the first-order autonomous differential system d dt ZðtÞ¼GðZðtÞÞ with Z ¼ðQ T ; _ Q T Þ T , the stability of the periodic solutions (70) can be deduced from the eigenvalues of the monodromy matrix associated with the fundamental matrix solution of the 2p-periodic variational linear differential system (see Ref. [11, p. 119 The monodromy matrix is computed over one period, using the four canonical basis vectors as initial conditions successively.The computations show that the periodic orbit (approximated by v Ã 1 ) associated with the first NM is stable on its invariant manifold and unstable in the phase space (two complex conjugate multipliers are outside the unit circle), whereas the periodic orbit (approximated by v Ã 2 ) associated with the  second NM is stable on its invariant manifold, as well as in the phase space (one multiplier lies on the unit circle and all the others are located inside the unit circle).The local and global stability are illustrated in Fig. 9 by solving numerically Eqs.(68).Motions are plotted with initial conditions (near the periodic orbit) on the invariant manifold, Fig. 9(a), and with initial conditions (near the periodic orbit) in the phase space outside the invariant manifold, Fig. 9(b).The two motions are attracted by the same periodic orbit approximated by v Ã 2 .

Conclusion
An amplitude-phase transformation procedure is described here for characterizing NMs (damped motion, periodic orbits, limit cycles, etc.) in the framework of invariant manifold theory.This system of formulation is suitable for use with a Galerkin-based computational procedure.Bifurcation analysis can be performed, and the existence and stability of periodic orbits on the associated invariant manifold can be studied from the two first-order coupled differential equations governing the amplitude and phase variables which describe the dynamics.If we look at the (generalized) damping rate function, it can be concluded that in some cases, periodic solutions may exist for the amplitude equation, reflecting the presence of stable or unstable ''limit cycles'' on the invariant manifold.Examples given here involve gyroscopic and/or damped nonlinear system equations and systems with auto-oscillations.It is worth noting that the present approach differs significantly from that described in Ref. [8] where, with an aim of improving the numerical procedure, an amplitude-phase transformation based on the (constant) frequency of the linearized system is used.

Fig. 2 .
Fig.2.Nonlinear mode which coincides for small values of v with the second mode of the underlying linear system of the non-gyroscopic conservative system: (a) X 1 versus v and f, (b) X 2 versus v and f, (c) O versus v and f and (d) period versus a.

Fig. 3 .
Fig. 3. Nonlinear mode which coincides for small values of v with the second mode of the underlying linear system of the gyroscopic conservative system: (a) X 1 versus v and f, (b) X 2 versus v and f, (c) O versus v and f and (d) period versus a.

Fig. 4 .
Fig.4.Nonlinear mode which coincides for small values of v with the second mode of the underlying linear system of the dissipative system with proportional damping: (a) X 1 versus v and f, (b) X 2 versus v and f, (c) O versus v and f and (d) x versus v and f.

c 11 ¼Fig. 5 .
Fig.5.Nonlinear mode which coincides for small values of v with the second mode of the underlying linear system of the dissipative system with non-proportional damping: (a) X 1 versus v and f, (b) X 2 versus v and f, (c) O versus v and f and (d) x versus v and f.

Fig. 9 .
Fig. 9. Illustration of the (a) local and (b) global stability of the limit cycle associated with the second nonlinear mode: limit cycle given by v Ã 2 (bold continuous lines), motions obtained by solving Eq. (68) with initial conditions inside the invariant manifold (a) (continuous line) and outside the invariant manifold (b) (continuous line).
and x are two scalar functions of the variables v and f, the frequency function and the damping function, respectively.As previously, we look for O and x in the form of even p-periodic functions with respect to f with Fourier series

Table 1
Period of the limit cycle v Ã 1 ¼ 1 for various values of e