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

Sergio Bellizzi, Robert Bouc

► To cite this version:

Sergio Bellizzi, Robert Bouc. An amplitude phase formulation for nonlinear modes and limit cycles analysis through invariant manifolds. Journal of Sound and Vibration, 2007, 300 (3-5), pp.896-915. 10.1016/j.jsv.2006.09.004. hal-00108694

HAL Id: hal-00108694

https://hal.science/hal-00108694v1

Submitted on 28 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An amplitude-phase formulation for nonlinear modes and limit cycles through invariant manifolds

S. Bellizzi*, R. Bouc

Laboratoire de Mécanique et d'Acoustique, CNRS, 31 Chemin Joseph Aiguier, F-13402 Marseille, France

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 defined in terms of amplitude, phase, frequency, damping coefficient 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 first-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.

## 1. Introduction

Nonlinear modes (NMs) are efficient tools for analyzing the behavior of vibrating mechanical systems [1–6]. In a series of papers by Shaw and Pierre [5–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

<footnote>*Corresponding author.</footnote>

<footnote>*E-mail address:* bellizzi@lma.cnrs-mrs.fr (S. Bellizzi).</footnote>
linear modal frequency in Ref. [8], takes computational advantage of the 2π-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 2π-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-of-freedom (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.

## 2. sdof Case

### 2.1. Undamped case with inertial nonlinearity

Let us consider an undamped sdof oscillator with inertial nonlinearity and a nonlinear restoring force. Let

$$T = \frac{1}{2}m(q)\dot{q}^2 \quad \text{and} \quad 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) = (\partial g/\partial q)(q)$, the Euler–Lagrange equation reads

$$\frac{d}{dt}(m(q)\dot{q}) - \frac{1}{2}m_q(q)\dot{q}^2 + U_q(q) = 0 \quad (1)$$

giving the equation of motion

$$m(q)\ddot{q} + \frac{1}{2}m_q(q)\dot{q}^2 + f(q) = 0, \quad (2)$$

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) \ge m_0 > 0$, $f(0) = 0$ and $qf(q) \ge 0$.

We will now establish that a periodic orbit of Eq. (2) can be expressed, either

$$\begin{aligned} q(t) &= a(\cos \phi(t) + b(a)), \\ \dot{q}(t) &= -a\Omega(a, \phi(t)) \sin \phi(t) \end{aligned} \quad (3)$$
with

$$ \dot{\phi}(t) = \Omega(a, \phi(t)), \quad \phi(0) = \phi \qquad (4) $$

or

$$ \begin{aligned} q(t) &= a(\cos \phi(t) + b(a, \phi(t))), \\ \dot{q}(t) &= a\Omega(a, \phi(t))(-\sin \phi(t) + b_{\phi}(a, \phi(t))) \end{aligned} \qquad (5) $$

with

$$ \dot{\phi}(t) = \Omega(a, \phi(t)), \quad \phi(0) = \phi. \qquad (6) $$

Numbers $a$ and $\phi$ 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 $\Omega$ and the bias function $b$ need to be found. These expressions differ from the classical polar transformation, since they involve a non-constant frequency $\Omega(a, \phi)$.

In expressions (3), $b$ is a scalar function which depends only on $a$, and $\Omega$ is a positive scalar function which depends on $a$ and $\phi$. As we shall see, $\Omega$ is an even $2\pi$-periodic continuous function with respect to $\phi$ and hence with Fourier series

$$ \Omega(a, \phi) = \sum_{n=0}^{\infty} \Omega_n(a) \cos n\phi. \qquad (7) $$

In expressions (5), the bias function $b$ depending on $a$ and $\phi$ will be sought for, arbitrarily, along with the frequency function $\Omega$, in the form of an even $\pi$-periodic function with respect to $\phi$ with Fourier series

$$ \begin{aligned} \Omega(a, \phi) &= \sum_{n=0}^{\infty} \Omega_{2n}(a) \cos 2n\phi, \\ b(a, \phi) &= \sum_{n=0}^{\infty} b_{2n}(a) \cos 2n\phi. \end{aligned} \qquad (8) $$

In both cases, it can be seen from Eqs. (4) and (6) that the period of the motion is given by

$$ T(a) = \int_{0}^{2\pi} \frac{1}{\Omega(a, \phi)} d\phi. \qquad (9) $$

After some algebraic and differential manipulations, substituting Eq. (3) into Eq. (2) gives

$$ \frac{d}{d\phi}(m(ax(\phi))\Omega^2 \sin^2 \phi) - \frac{2}{a} \sin \phi f(ax(\phi)) = 0, \qquad (10) $$

where $x(\phi) = \cos \phi + b(a)$ yielding

$$ \Omega^2(a, \phi) = \frac{2 \int_0^\phi \sin \sigma f(ax(\sigma)) d\sigma}{am(ax(\phi)) \sin^2 \phi} = \frac{2(U(a + ab(a)) - U(a \cos \phi + ab(a)))}{a^2 m(ax(\phi)) \sin^2 \phi}, \qquad (11) $$

with the constraint equation

$$ \int_{0}^{\pi} \sin \sigma f(ax(\sigma)) d\sigma = 0, \qquad (12) $$

i.e.

$$ U(a(1 + b(a))) - U(a(-1 + b(a))) = 0. \qquad (13) $$

Eq. (13) defines the bias function $b(a)$ and, together with the assumption $qf(q) \ge 0$, means that $\Omega^2$, given by Eq. (11), is a positive continuous $2\pi$-periodic function with respect to $\phi$ 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

$$ \frac{d}{d\phi}(m(ax(\phi))\Omega^2 \sin^2 \phi) - \frac{2}{a}\sin \phi f(ax(\phi)) \\ - \sin \phi \left( \frac{d}{d\phi}(m(ax(\phi))\Omega^2 b_{\phi}) + m(ax(\phi))\Omega^2 b_{\phi\phi} \right) = 0, \quad (14) $$

where now $x(\phi) = \cos \phi + b(a, \phi)$. 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 $\Omega^2$ and $b$. Nevertheless we look for even $\pi$-periodic functions (see Eq. (8)). Consequently, using the decomposition

$$ \begin{aligned} m(ax(\phi)) &= \sum_{n=0}^{\infty} m_{2n}(a) \cos 2n\phi + \sum_{n=0}^{\infty} m_{2n+1}(a) \cos(2n+1)\phi \\ &= m^e(\phi) + m^o(\phi), \end{aligned} $$

$$ \begin{aligned} f(ax(\phi)) &= \sum_{n=0}^{\infty} f_{2n}(a) \cos 2n\phi + \sum_{n=0}^{\infty} f_{2n+1}(a) \cos(2n+1)\phi \\ &= f^e(\phi) + f^o(\phi), \end{aligned} $$

where $(\cdot)^e$ (respectively $(\cdot)^o$) denotes the even (respectively, odd) cosine terms in Fourier series, Eq. (14) can be re-written as

$$ \frac{d}{d\phi}(m^e(\phi)\Omega^2 \sin^2 \phi) - \frac{2}{a}\sin \phi f^o(\phi) \\ - \sin \phi \left( \frac{d}{d\phi}(m^o(\phi)\Omega^2 b_{\phi}) + m^o(\phi)\Omega^2 b_{\phi\phi} \right) = 0, \quad (15) $$

$$ \frac{d}{d\phi}(m^o(\phi)\Omega^2 \sin^2 \phi) - \frac{2}{a}\sin \phi f^e(\phi) \\ - \sin \phi \left( \frac{d}{d\phi}(m^e(\phi)\Omega^2 b_{\phi}) + m^e(\phi)\Omega^2 b_{\phi\phi} \right) = 0, \quad (16) $$

where Eq. (15) (respectively, Eq. (16)) stands for the odd (respectively, even) sine terms of Eq. (14).

If $m^o \equiv 0$ and $f^e \equiv 0$, then $b(a, \phi) = 0$ is a solution of Eq. (16), and Eq. (15) reduces to Eq. (10), yielding the same expression for Eqs. (3) and (5) and $\Omega$ is $\pi$-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.

## 2.2. Odd restoring force with linear damping

Consider a sdof oscillator with linear damping, unit mass and a nonlinear restoring force

$$ \ddot{q} + c\dot{q} + f(q) = 0. \quad (17) $$

We will assume that the function $f$ is regular enough and, for the sake of simplicity, that $f(q) = -f(-q)$.

If $c \neq 0$, 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:

$$ q(t) = v(t) \cos \phi(t) \quad (18) $$

with

$$ \begin{aligned} \dot{\phi}(t) &= \Omega(v(t), \phi(t)), \\ \dot{v}(t) &= v(t)\xi(v(t), \phi(t)), \end{aligned} \quad \phi(0) = \varphi, \quad v(0) = a, \quad (19) $$
where $\Omega$ and $\xi$ are two scalar functions of the variables $v$ and $\phi$, the frequency function and the damping function, respectively. As previously, we look for $\Omega$ and $\xi$ in the form of even $\pi$-periodic functions with respect to $\phi$ with Fourier series

$$ 
\begin{aligned}
\Omega(v, \phi) &= \sum_{n=0}^{\infty} \Omega_{2n}(v) \cos 2n\phi, \\
\xi(v, \phi) &= \sum_{n=0}^{\infty} \xi_{2n}(v) \cos 2n\phi.
\end{aligned} 
\qquad (20) 
$$

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 $\phi$ for the two unknown functions $\Omega$ and $\xi$

$$ \frac{\partial \xi}{\partial \phi} \cos \phi - 2\xi \sin \phi - \xi v \frac{1}{2\Omega^2} \frac{\partial \Omega^2}{\partial v} \sin \phi = c \sin \phi, \qquad (21) $$

$$ \frac{1}{2} \frac{\partial \Omega^2}{\partial \phi} \sin \phi + \Omega^2 \cos \phi - \frac{1}{v} f(v \cos \phi) = \cos \phi \left( \xi^2 + \frac{1}{2} v \frac{\partial \xi^2}{\partial v} + c \xi \right). \qquad (22) $$

If $f(q) = kq$, then $\Omega^2 = k - (c^2/4)$ and $\xi = -c/2$ as was to be expected. Otherwise, it can be established that in the case of positive (respectively, negative) damping, $c>0$ (respectively, $c<0$), we have $\xi \le 0$ (respectively, $\xi \ge 0$) 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.

### 2.3. Application to limit cycles

From Eq. (19), it follows that

$$ \frac{dv}{d\phi} = v\tau(v, \phi), \qquad (23) $$

where $\tau(v, \phi) = (\xi(v, \phi)/\Omega(v, \phi))$ can be viewed as a “generalized damping rate function”.

Since the right-hand side of Eq. (23) is $\pi$-periodic with respect to $\phi$, periodic solutions

$$ v^*(\phi) = v^*(\phi + \pi) $$

may exist with some $\xi$ and $\Omega$ (one necessary condition being that $\tau(v, \phi)$ does not keep a constant sign). It follows that

$$ q(t) = v^*(\phi(t)) \cos \phi(t) $$

will be a $T$-periodic function with period

$$ T = \int_{0}^{2\pi} \frac{d\phi}{\Omega(v^*(\phi), \phi)}. \qquad (24) $$

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

$$ \ddot{q} - \varepsilon(1 - \frac{22}{3}q^2 + 8q^4 - \frac{32}{15}q^6)\dot{q} + q = 0, \qquad (25) $$

where $\varepsilon$ is a positive (small) parameter.
Substituting Eqs. (18) and (19) into Eq. (25) gives similar equations to Eqs. (21) and (22)

$$ \frac{\partial \xi}{\partial \phi} \cos \phi - 2\xi \sin \phi - \xi v \frac{1}{2\Omega^2} \frac{\partial \Omega^2}{\partial v} \sin \phi = \epsilon p(v \cos \phi) \sin \phi, \quad (26) $$

$$ \frac{1}{2} \frac{\partial \Omega^2}{\partial \phi} \sin \phi + \Omega^2 \cos \phi - \cos \phi = \cos \phi \left( \xi^2 + \frac{1}{2} v \frac{\partial \xi^2}{\partial v} + \epsilon p(v \cos \phi) \xi \right), \quad (27) $$

where $p(q) = -(1 - \frac{22}{3}q^2 + 8q^4 - \frac{32}{15}q^6)$. One can check that solutions of Eqs. (26) and (27) are such that

$$ \xi = O(\epsilon), \quad \Omega^2 = 1 + O(\epsilon^2) \quad \text{and} \quad \xi \frac{1}{\Omega^2} \frac{\partial \Omega^2}{\partial v} = O(\epsilon^3). $$

It follows from Eq. (26), that

$$ \xi(v, \phi) = \xi_0(v) + \xi_2(v) \cos 2\phi + \xi_4(v) \cos 4\phi + \xi_6 \cos 6\phi + O(\epsilon^3) \quad (28) $$

with

$$ \xi_0(v) = -\frac{\epsilon}{12}(-6 + 11v^2 - 6v^4 + v^6), \quad \xi_2(v) = -\frac{\epsilon}{24}(22v^2 - 16v^4 + 3v^6), $$

$$ \xi_4(v) = -\frac{\epsilon}{60}(-10v^4 + 3v^6), \quad \xi_6(v) = -\frac{\epsilon}{120}v^6 $$

and, from Eq. (27), we deduce that

$$ \Omega^2(v, \phi) = \sum_{n=0}^{6} \Omega_{2n}^2(v) \cos 2n\phi + O(\epsilon^4) \quad (29) $$

with

$$ \Omega_0^2(v) = 1 + \epsilon^2 \left( -\frac{1}{4} + \frac{11}{4}v^2 - \frac{845}{96}v^4 + \frac{203}{18}v^6 - \frac{217}{32}v^8 + \frac{77}{40}v^{10} - \frac{1001}{4800}v^{12} \right) $$

and, for $n > 0$, $\Omega_{2n}^2 = O(\epsilon^2)$.

Substituting Eqs. (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 $\xi = \epsilon\tilde{\xi}$ and hence $\tau = \epsilon\tilde{\tau}$)

$$ \frac{dv}{d\phi} = v(\tau(v)), \quad \text{where} \quad \langle \tau(\cdot) \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} \frac{\xi(\cdot, \phi)}{\Omega(\cdot, \phi)} d\phi, \quad (30) $$

yielding

$$ \langle \tau(v) \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} \xi(v, \phi) d\phi + O(\epsilon^2) = \xi_0(v) + O(\epsilon^2). $$

The amplitudes of the limit cycles are given by the three positive roots of the equation $\xi_0(v) = 0$ up to $O(\epsilon^2)$. The first root $v_1^* = 1$ gives a stable limit cycle ($\partial\xi_0(1)/\partial v < 0$), the second root $v_2^* = \sqrt{2}$ gives an unstable limit cycle ($\partial\xi_0(\sqrt{2})/\partial v > 0$), and the third root $v_3^* = \sqrt{3}$ gives a stable limit cycle ($\partial\xi_0(\sqrt{3})/\partial v < 0$). 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

$$ T_i = \frac{2\pi}{\Omega_0(v_i^*)} + O(\epsilon^2). \quad (31) $$

Table 1 compares, for various values of the parameter $\epsilon$, 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.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Phase portrait $\epsilon = 0.5$, $q(0) = \sqrt{2} - 0.1$ (continuous line) and $q(0) = \sqrt{2} + 0.1$ (dashed line).</figcaption></figure>

<table><caption>Table 1
Period of the limit cycle $v_1^* = 1$ for various values of $\epsilon$</caption><thead><tr><th></th><th>ε = 0.01</th><th>ε = 0.1</th><th>ε = 0.5</th><th>ε = 0.8</th><th>ε = 1</th></tr></thead><tbody><tr><td>T given by Eq. (31)</td><td>6.2832</td><td>6.2860</td><td>6.3543</td><td>6.4710</td><td>6.5833</td></tr><tr><td>T given by a shooting method</td><td>6.2832</td><td>6.2849</td><td>6.3263</td><td>6.3938</td><td>6.4564</td></tr></tbody></table>

### 3. mdof Systems

We now consider a mdof system of the form

$$[\mathbf{M}]\ddot{\mathbf{Q}}(t) + \mathbf{F}(\dot{\mathbf{Q}}(t), \mathbf{Q}(t)) = \mathbf{0}, \quad (32)$$

where $\mathbf{Q}(t)$ is an $n$-vector, $[\mathbf{M}]$ is a non-singular symmetric $n \times n$-matrix and $\mathbf{F}$ is a (sufficiently regular) vector function with dimension $n$ such that $\mathbf{F}(\mathbf{0}, \mathbf{0}) = \mathbf{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.

#### 3.1. Undamped case

In a previous work [9] we considered the case in which the nonlinear force function $\mathbf{F}$ does not depend on $\dot{\mathbf{Q}}$. We briefly summarize here the main results.

Assume firstly that $\mathbf{F}$ is an odd function ($\mathbf{F}(\mathbf{Q}) = -\mathbf{F}(-\mathbf{Q})$). A normal mode is defined from a family of periodic solutions to Eq. (32) sought in the form

$$\mathbf{Q}(t) = a\mathbf{\Psi}(a, \phi(t)) \cos \phi(t), \quad (33)$$

where the scalar function $\phi$ satisfies the differential equation

$$\dot{\phi}(t) = \Omega(a, \phi(t)), \quad \phi(0) = \phi. \quad (34)$$

By analogy with the linear case (where Eqs. (33) and (34) hold with constant $\mathbf{\Psi}$ and $\Omega$), the vector function, $\mathbf{\Psi}$ and the scalar function $\Omega$, will be referred to here as the modal vector and the resonance frequency of the nonlinear normal mode, respectively, and the scalar function $\phi$, will be referred to here as the total phase of the motion. The one-dimensional differential equation (34) governing the total phase motion $\phi$, will define the dynamics of the periodic response and the scalar quantities $a$ ($a > 0$) and $\phi$ ($\phi \in [0, 2\pi]$) will set the initial conditions of the vibration in the phase space. To ensure that the parameter $a$ appropriately characterizes the
amplitude of the vibration in mode motion, a normalization condition on the modal vector $\Psi$ is required. One can set, for example, for all $(a, \phi)$,

$$ \Psi^T(a, \phi)[\mathbf{M}]\Psi(a, \phi) = 1. \quad (35) $$

Substituting Eq. (33) into Eq. (32) and eliminating the time variable in the resulting equation by using Eq. (34), Eq. (32) reduces to the following differential equation in the variable $\phi$:

$$ [\mathbf{M}]\Psi\left(\Omega^2 \cos \phi + \frac{1}{2}(\Omega^2)_{\phi} \sin \phi\right) = \mathbf{L}(\Omega^2, \Psi; \phi) + \frac{1}{a}\mathbf{F}(\Psi a \cos \phi), \quad (36) $$

where

$$ \mathbf{L}(\Omega^2, \Psi; \phi) = \Omega^2 \cos \phi [\mathbf{M}] \Psi_{\phi\phi} - 2\Omega^2 \sin \phi [\mathbf{M}] \Psi_{\phi} + \frac{1}{2}(\Omega^2)_{\phi} \cos \phi [\mathbf{M}] \Psi_{\phi} \quad (37) $$

and $(\cdot)_{\phi}$ denotes partial differentiation with respect to $\phi$. The differential rule $(\Omega^2)_{\phi} = 2\Omega\Omega_{\phi}$ has been used to work with the unknown function $\Omega^2$ in place of $\Omega$.

Eq. (36), which is 2$\pi$-periodic with respect to $\phi$, together with the normalization condition (35), can be viewed as a nonlinear 2$\pi$-periodic eigenvalue–eigenvector problem. Thus, the modal vector $\Psi$ and the frequency $\Omega$ will be searched for as 2$\pi$-periodic functions with respect to $\phi$ for fixed $a$. More specifically, due to the odd symmetry of the function $F$, and the quadratic form (in $\Psi$) of the normalization constraint (35), these functions will be searched for as even periodic functions with respect to $\phi$ with period $\pi$ and hence with Fourier series containing only cosine terms of even order. It follows that the Fourier series of $Q(t)$ in terms of $\phi(t)$ will contain only cosines terms of odd order. It can be established that if $Q^T F(Q) \ge 0$ (which is weaker than the monotony assumption made in Ref. [9]) necessarily leads to $\Omega^2 \ge 0$. The period of the vibration in mode motion depends only on the amplitude $a$ and follows easily from Eq. (34):

$$ T(a) = \int_{0}^{2\pi} \frac{1}{\Omega(a, \phi)} d\phi. \quad (38) $$

Eq. (33) together with

$$ \dot{\mathbf{Q}}(t) = a\Omega(a, \phi(t))(\Psi_{\phi}(a, \phi(t))\cos\phi(t)) - \Psi(a, \phi(t))\sin\phi(t), \quad (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 $\phi$ to be independent variables.

If $F(\mathbf{Q}) \neq -F(-\mathbf{Q})$, due to the presence of quadratic terms for example, Eq.(33) must be replaced by [9]

$$ \mathbf{Q}(t) = a(\Psi(a, \phi(t))\cos(\phi(t)) + \mathbf{B}(a, \phi(t))), \quad (40) $$

where the term $B(a, \phi)$, such that $B(a, \phi) = B(a, -\phi) = B(a, \phi + \pi)$, 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 $(\mathbf{Q}, \dot{\mathbf{Q}})$ in the form

$$ \begin{aligned} \mathbf{Q}(t) &= a\mathbf{X}(a, \phi(t)), \\ \dot{\mathbf{Q}}(t) &= a\mathbf{Y}(a, \phi(t)), \end{aligned} \quad (41) $$

where $\phi(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 $\phi$ for the unknowns $\mathbf{X}, \mathbf{Y}, \Omega$,

$$ \mathbf{X}_{\phi}\Omega = \mathbf{Y}, \quad (42) $$

$$ [\mathbf{M}]\mathbf{Y}_{\phi}\Omega + \frac{1}{a}\mathbf{F}(a\mathbf{X}) = \mathbf{0} \quad (43) $$
for which 2π-periodic solutions are planned. To complete the above equations a normalization condition must be added. We use the following result: any regular, even, 2π-periodic vector function X(φ) can always be written X(φ) = X<sup>oc</sup>(φ) + X<sup>ec</sup>(φ), where X<sup>oc</sup> (resp., X<sup>ec</sup>) stands for the odd cosine (resp., even cosine) terms in the Fourier series of X(φ). Moreover, there exists one and only one even and π-periodic vector function in φ, Ψ<sub>c</sub>(φ), such that X<sup>oc</sup>(φ) = Ψ<sub>c</sub>(φ) cos φ. It follows that the normalization constraint (35) can again be retained and rewritten as

$$
X^{\text{oc}T}[M]X^{\text{oc}} = \cos^2 \phi. \tag{44}
$$

Finally, we seek periodic solutions of Eqs. (42)–(44) such that

$$
X(a, \phi) = X(a, \phi + 2\pi) = X(a, -\phi), \quad (45)
$$

$$
\mathbf{Y}(a, \phi) = \mathbf{Y}(a, \phi + 2\pi) = -\mathbf{Y}(a, -\phi), \qquad (46)
$$

$$
\Omega(a, \phi) = \Omega(a, -\phi) = \Omega(a, \phi + \pi). \tag{47}
$$

3.2. 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),dotQ(t)), are related to a
single pair of amplitude and phase variables, (v(t), φ(t)), according to

$$
\begin{align}
\mathbf{Q}(t) &= v(t)\mathbf{X}(v(t), \phi(t)), \notag \\
\dot{\mathbf{Q}}(t) &= v(t)\mathbf{Y}(v(t), \phi(t)), \tag{48}
\end{align}
$$

where X and Y are *n*-vector functions. The amplitude and phase variables are governed by the two first-order differential equations

$$
\begin{align}
\dot{\phi}(t) &= \Omega(v(t), \phi(t)), \nonumber \\
\dot{v}(t) &= v(t) \xi(v(t), \phi(t)), \tag{49}
\end{align}
\[
\phi(0) = \varphi, \quad v(0) = a,
\]
$$

where Ω and ξ are scalar functions. Numbers φ (∈ [0, 2π]) and a (>0) are two given constants which set the initial conditions of the motion.

As in Section 3.1, the four functions **X**, **Y**, Ω and ξ will define a modal motion for Eq. (32), *provided that they*
*are 2π-periodic in the phase variable* *φ*. Again, the two-dimensional invariant manifold associated with the NM
*in the phase space reads*

$$
\begin{align}
\mathbf{Q} &= a\mathbf{X}(a, \varphi), \notag \\
\dot{\mathbf{Q}} &= a\mathbf{Y}(a, \varphi), \tag{50}
\end{align}
$$

while the two first-order differential equations (49) govern the dynamics of the motion on the invariant manifold.

The two scalar functions, Ω (the frequency function) and ξ (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, φ),

$$
(X + vX_v)\xi + X_\phi\Omega = Y, \tag{51}
$$

$$
[\mathbf{M}](\mathbf{Y} + v\mathbf{Y}_v)\xi + [\mathbf{M}]\mathbf{Y}_\phi\Omega + \frac{1}{v}\mathbf{F}(v\mathbf{Y}, v\mathbf{X}) = \mathbf{0}. \quad (52)
$$
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 **φ** and the two unknown scalar functions **Ω** and **ξ**.

To characterize the four unknown functions **X**, **Y**, **Ω** and **ξ**, it is necessary to add two constraint scalar equations to Eqs. (51) and (52) (normalization conditions). Due to the 2π-periodicity in the phase variable **φ**, the function **X** can always be broken up according to

$$ \mathbf{X} = \mathbf{X}^{\mathrm{oc}} + \mathbf{X}^{\mathrm{ec}} + \mathbf{X}^{\mathrm{os}} + \mathbf{X}^{\mathrm{es}}, $$

where (.)<sup>oc</sup> ((.)<sup>ec</sup>, (.)<sup>os</sup>, (.)<sup>es</sup>, 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 π-periodic vector functions **Ψ**<sub>c</sub>, **Ψ**<sub>s</sub> such that

$$ \mathbf{X}^{\mathrm{oc}}(v, \phi) = \mathbf{\Psi}_{\mathrm{c}}(v, \phi) \cos \phi, \quad \mathbf{X}^{\mathrm{os}}(v, \phi) = \mathbf{\Psi}_{\mathrm{s}}(v, \phi) \sin \phi. $$

It follows that we can adopt the normalization condition (usually set for linear systems):

$$ \mathbf{\Psi}_{\mathrm{c}}^{\mathrm{T}}[\mathbf{M}]\mathbf{\Psi}_{\mathrm{c}} + \mathbf{\Psi}_{\mathrm{s}}^{\mathrm{T}}[\mathbf{M}]\mathbf{\Psi}_{\mathrm{s}} = 1, $$

$$ \mathbf{\Psi}_{\mathrm{c}}^{\mathrm{T}}[\mathbf{M}]\mathbf{\Psi}_{\mathrm{s}} = 0, $$

which can be rewritten as

$$ \sin^2 \phi \mathbf{X}^{\mathrm{ocT}}[\mathbf{M}]\mathbf{X}^{\mathrm{oc}} + \cos^2 \phi \mathbf{X}^{\mathrm{osT}}[\mathbf{M}]\mathbf{X}^{\mathrm{os}} = \sin^2 \phi \cos^2 \phi, \quad (53) $$

$$ \mathbf{X}^{\mathrm{ocT}}[\mathbf{M}]\mathbf{X}^{\mathrm{osT}} \sin \phi \cos \phi = 0. \quad (54) $$

Finally, a NM of the system (32) is obtained by solving Eqs. (51)–(54) for the four functions **X**, **Y**, **Ω** and **ξ**, recalling that

$$ \mathbf{X}(v, \phi) = \mathbf{X}(v, \phi + 2\pi), \quad \mathbf{Y}(v, \phi) = \mathbf{Y}(v, \phi + 2\pi), \quad (55) $$

$$ \Omega(v, \phi) = \Omega(v, -\phi) = \Omega(v, \phi + \pi), \quad \xi(v, \phi) = \xi(v, -\phi) = \xi(v, \phi + \pi). \quad (56) $$

If **F**(Q̇, Q) = −**F**(-Q̇, −Q) then v**X**(v, **φ**) takes the form

$$ v\mathbf{X}(v, \phi) = v(\mathbf{\Psi}_{\mathrm{c}}(v, \phi) \cos \phi - \mathbf{\Psi}_{\mathrm{s}}(v, \phi) \sin \phi), $$

which can be interpreted in terms of amplitude modulation coefficient (v), eigenvectors (Ψ<sub>c</sub>, Ψ<sub>s</sub>) and frequency (Ω), like for a linear system.

## 4. Computing the nonlinear modes

The challenge is now how to solve the PDEs (51)–(54) for the four unknown functions **X**, **Y**, **Ω** and **ξ**. These equations have been formally established for (v, φ) ∈ ℝ<sup>+</sup> × [0, 2π]. In the approximation procedure, the domain over which the PDEs will be solved will be limited to [0, v<sub>max</sub>] × [0, 2π], where v<sub>max</sub> is a positive constant.

In the case of small oscillation responses (i.e. |v| ≪ 1), denoting by [∂**F**<sub>Q</sub>(., .)] and by [∂**F**<sub>Q</sub>(., .)] the Jacobian matrix of **F** with respect to Q̇ and Q, the nonlinear system (32) can be replaced by the linearized one

$$ [\mathbf{M}]\ddot{\mathbf{Q}}(t) + [\partial\mathbf{F}_{\dot{\mathbf{Q}}}(\mathbf{0}, \mathbf{0})]\dot{\mathbf{Q}} + [\partial\mathbf{F}_{\mathbf{Q}}(\mathbf{0}, \mathbf{0})]\mathbf{Q} = \mathbf{0}, \quad (57) $$

which gives *n* linear modes. These linear modes solve the PDEs (51)–(54) by introducing into Eq. (52) the following approximation:

$$ \frac{\mathbf{F}(v\mathbf{X}(v, \phi), v\mathbf{Y}(v, \phi))}{v} \approx [\partial\mathbf{F}_{\dot{\mathbf{Q}}}(\mathbf{0}, \mathbf{0})]\mathbf{Y}(v, \phi) + [\partial\mathbf{F}_{\mathbf{Q}}(\mathbf{0}, \mathbf{0})]\mathbf{X}(v, \phi). $$
Each one of these linear modes can be written

$$
\begin{align}
\mathbf{X}(v, \phi) &= \mathbf{\Psi}_c \cos \phi - \mathbf{\Psi}_s \sin \phi, \nonumber \\
\mathbf{Y}(v, \phi) &= (\xi \mathbf{\Psi}_c - \Omega \mathbf{\Psi}_s) \cos \phi - (\xi \mathbf{\Psi}_s + \Omega \mathbf{\Psi}_c) \sin \phi, \tag{58} \\
\Omega(v, \phi) &= \omega \quad \text{and} \quad \xi(v, \phi) = \eta. \nonumber
\end{align}
$$

where the complex frequency $\lambda = \eta + i\omega$ ($i = \sqrt{-1}$), solves

$$
\det([\mathbf{M}]\lambda^2 + [\partial\mathbf{F}_{\dot{Q}}(\mathbf{0}, \mathbf{0})]\lambda + [\partial\mathbf{F}_{Q}(\mathbf{0}, \mathbf{0})]) = 0. \quad (59)
$$

The eigenvector, $\mathbf{\Psi} = \mathbf{\Psi}_c + i\mathbf{\Psi}_s$, is next given by solving the two following equations:

$$
([\mathbf{M}](\eta^2 - \omega^2) + \eta[\partial\mathbf{F}_Q(\mathbf{0}, \mathbf{0})] + [\partial\mathbf{F}_Q(\mathbf{0}, \mathbf{0})])\mathbf{\Psi}_c - (2\omega\eta[\mathbf{M}] + \omega[\partial\mathbf{F}_Q(\mathbf{0}, \mathbf{0})])\mathbf{\Psi}_s = 0, \quad (60)
$$

$$
(2\omega\eta[\mathbf{M}] + \omega[\partial\mathbf{F}_{\dot{Q}}(\mathbf{0}, \mathbf{0})])\mathbf{\Psi}_c + ([\mathbf{M}](\eta^2 - \omega^2) + \eta[\partial\mathbf{F}_{\dot{Q}}(\mathbf{0}, \mathbf{0})] + [\partial\mathbf{F}_{\dot{Q}}(\mathbf{0}, \mathbf{0})])\mathbf{\Psi}_s = 0 \quad (61)
$$

under the normalization constraints

$$
\mathbf{\Psi}_c^T [\mathbf{M}] \mathbf{\Psi}_c + \mathbf{\Psi}_s^T [\mathbf{M}] \mathbf{\Psi}_s = 1, \quad (62)
$$

$$
\mathbf{\Psi}_c^T [\mathbf{M}] \mathbf{\Psi}_s = 0. \tag{63}
$$

We focus now on calculating the NMs which coincide at small values of $v$ (i.e. $v \approx 0$) with one of the modes
of the linear system (57).

Since the unknown functions are 2π-periodic and π-periodic with respect to the variable φ (see Eqs. (55) and (56)), a Galerkin discretization procedure based on the truncated Fourier series with respect to the variable φ 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

$$
\begin{align}
\mathbf{X}^{\text{oc}}(v, \phi) &\approx \mathbf{X}_{\text{lin}}^{\text{c}}(\phi) + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{X}_{p, 2k+1}^{\text{c}} v^p \cos(2k+1)\phi, \nonumber \\
\mathbf{X}^{\text{ec}}(v, \phi) &\approx \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{X}_{p, 2k}^{\text{c}} v^p \cos 2k\phi, \nonumber \\
\mathbf{X}^{\text{os}}(v, \phi) &\approx \mathbf{X}_{\text{lin}}^{\text{s}}(\phi) + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{X}_{p, 2k+1}^{\text{s}} v^p \sin(2k+1)\phi, \nonumber \\
\mathbf{X}^{\text{es}}(v, \phi) &\approx \sum_{p=1}^{N_v} \sum_{k=1}^{N_{\phi}} \mathbf{X}_{p, 2k}^{\text{s}} v^p \sin 2k\phi, \nonumber
\end{align}
$$

$$
\begin{align}
\mathbf{Y}^{\text{oc}}(v, \phi) &\approx \mathbf{Y}_{\text{lin}}^{\text{c}}(\phi) + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{Y}_{p, 2k+1}^{\text{c}} v^p \cos(2k+1)\phi, \nonumber \\
\mathbf{Y}^{\text{ec}}(v, \phi) &\approx \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{Y}_{p, 2k}^{\text{c}} v^p \cos 2k\phi, \nonumber \\
\mathbf{Y}^{\text{os}}(v, \phi) &\approx \mathbf{Y}_{\text{lin}}^{\text{s}}(\phi) + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \mathbf{Y}_{p, 2k+1}^{\text{s}} v^p \sin(2k+1)\phi, \tag{65} \\
\mathbf{Y}^{\text{es}}(v, \phi) &\approx \sum_{p=1}^{N_v} \sum_{k=1}^{N_{\phi}} \mathbf{Y}_{p, 2k}^{\text{s}} v^p \sin 2k\phi, \nonumber
\end{align}
$$
$$ 
\begin{aligned}
\Omega(v, \phi) &\approx \Omega_{\text{lin}} + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \Omega_{p,2k} v^p \cos 2k\phi, \\
\xi(v, \phi) &\approx \xi_{\text{lin}} + \sum_{p=1}^{N_v} \sum_{k=0}^{N_{\phi}} \xi_{p,2k} v^p \cos 2k\phi,
\end{aligned}
\qquad (66) $$

where $N_v \ge 1$ and $N_{\phi} \ge 0$ denote the order of the truncated series and $X_{p,2k+1}^c, X_{p,2k}^c, \dots, \xi_{p,2k}$ denote the unknown expansion coefficients which have to be determined. The vector function $X_{\text{lin}}^c, X_{\text{lin}}^s, Y_{\text{lin}}^c$ and $Y_{\text{lin}}^s$ and the constants $\Omega_{\text{lin}}$ and $\xi_{\text{lin}}$ which set the linear mode under investigation are calculated from Eqs. (58)–(63).

Substituting Eqs. (64)–(66) into (51)–(54) and applying a Galerkin projection based on the scalar product

$$ (f,g) = \int_{0}^{v_{\text{max}}} \int_{0}^{2\pi} f(v, \phi)g(v, \phi) \, dv \, d\phi $$

and the basis functions $v^p \cos m\phi, v^p \sin m\phi$ (in Eqs. (51) and (52)) and $v^p \cos 2m\phi$ (in Eqs. (53) and (54)) leads to a set of $4N_v(2N_{\phi}+1)n + 2N_vN_{\phi}$ nonlinear algebraic equations in the unknown $X_{p,2k+1}^c, X_{p,2k}^c, \dots, \xi_{p,2k}$. For fixed $N_v$ and $N_{\phi}$, 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^{\text{ec}} \equiv 0, X^{\text{es}} \equiv 0, Y^{\text{ec}} \equiv 0, Y^{\text{es}} \equiv 0$.

# 5. Examples

## 5.1. Linear damping and nonlinear restoring forces

Let us consider the two-dimensional mechanical system presented in Refs. [8,9]

$$ 
\begin{aligned}
\ddot{q}_1 + c_{11}\dot{q}_1 + c_{12}\dot{q}_2 + \omega_1^2 q_1 + 0.405q_1^3 + 1.34q_1^2q_2 + 1.51q_1q_2^2 + 0.349q_2^3 &= 0, \\
\ddot{q}_2 + c_{21}\dot{q}_1 + c_{22}\dot{q}_2 + \omega_2^2 q_2 + 0.448q_1^3 + 1.51q_1^2q_2 + 1.05q_1q_2^2 + 4.580q_2^3 &= 0
\end{aligned}
\qquad (67) $$

with $\omega_1 = 0.689$ and $\omega_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_{\text{max}}] \times [0, 2\pi]$ with $v_{\text{max}} = 1.9$ as described in Section 4. Due to the symmetrical properties of the system, we have $X^{\text{ec}} \equiv 0, X^{\text{es}} \equiv 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_{\text{max}}]$ defined by

$$ \rho^2(T_{\text{max}}) = \frac{\int_0^{T_{\text{max}}} \|\mathbf{Q}^a(t) - \mathbf{Q}^e(t)\|^2 dt}{\int_0^{T_{\text{max}}} \|\mathbf{Q}^e(t)\|^2 dt}, $$

where $\mathbf{Q} = (q_1, q_2)$. $\mathbf{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_{\text{max}}$ and $\varphi = 0$. $\mathbf{Q}^e$ denotes the corresponding motion obtained by solving the equation of motion Eq. (67) with the initial condition $\mathbf{Q}^e(0) = \mathbf{Q}^a(0)$ and $\dot{\mathbf{Q}}^e(0) = \dot{\mathbf{Q}}^a(0)$. The range of integration $T_{\text{max}}$ depends on the behavior of the modal motion under investigation.

### 5.1.1. Non-gyroscopic conservative system

Parameter values: $c_{11} = c_{12} = c_{21} = c_{22} = 0$. Linear mode under consideration (see Section 4): $\eta = 0$, $\omega = 3.244$, $\Psi^c = (0, 0)^T$ and $\Psi^s = (0, 1)^T$.
<figure><img src="image_3.png" /><figcaption>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 $\phi$, (b) $X_2$ versus $v$ and $\phi$, (c) $\Omega$ versus $v$ and $\phi$ and (d) period versus $a$.</figcaption></figure>

The corresponding NM is shown in Figs. 2 (a)–(d). All the coefficients in the $\xi$-expansion (Eq. (66)) as well as those in the $X^{\text{oc}}$-expansion (Eq. (64)) were found to be equal to zero. An $X$-expansion with only cosine terms was obtained upon taking $\Psi^c = (0, 1)^T$ and $\Psi^s = (0, 0)^T$. This result is in agreement with the amplitude-phase approach performed in Ref. [9], where the functions $\xi$ and $X^{\text{os}}$ were both taken to be equal to zero. Since function $\xi$ 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_{\text{max}} = T(v_{\text{max}})$ (given by Eq. (38)) with different values of $N_v$ and $N_\phi$ are given in Table 2. Orders $N_v = 4$ and $N_\phi = 2$, give sufficiently accurate approximations over $[0, 2\pi] \times [0, v_{\text{max}}]$.

### 5.1.2. Gyroscopic conservative system

Parameter values: $c_{11} = c_{22} = 0$, $c_{12} = -c_{21} = 1/3$. Linear mode under consideration: $\eta = 0$, $\omega = 3.262$, $\Psi^c = (-0.106, 0)^T$ and $\Psi^s = (0, 0.994)^T$.

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
<table><caption>Table 2
Normalized mean square error</caption><thead><tr><th></th><th>N<sub>v</sub> = 2<br>N<sub>φ</sub> = 0</th><th>N<sub>v</sub> = 2<br>N<sub>φ</sub> = 1</th><th>N<sub>v</sub> = 2<br>N<sub>φ</sub> = 2</th><th>N<sub>v</sub> = 4<br>N<sub>φ</sub> = 0</th><th>N<sub>v</sub> = 4<br>N<sub>φ</sub> = 1</th><th>N<sub>v</sub> = 4<br>N<sub>φ</sub> = 2</th></tr></thead><tbody><tr><th>ρ</th><td>0.10277</td><td>0.04051</td><td>0.04183</td><td>0.13086</td><td>0.00651</td><td>0.00473</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>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 $\phi$, (b) $X_2$ versus $v$ and $\phi$, (c) $\Omega$ versus $v$ and $\phi$ and (d) period versus $a$.</figcaption></figure>

coefficients in the $\xi$-expansion (Eq. (66)) were not found to be both equal to zero, but to have small values. It was checked that imposing $\xi \equiv 0$ yields the same result. The normalized mean square error was $\rho(T(v_{max})) = 0.00433$ with $N_v = 4$ and $N_\phi = 2$.

### 5.1.3. Dissipative system with proportional damping

Parameter values: $c_{12} = c_{21} = 0$, $c_{ii} = 2\tau\omega_i$ with $\tau = 0.2$. Linear mode under consideration: $\eta = -0.649$, $\omega = 3.178$, $\Psi^c = (0, 0)^T$ and $\Psi^s = (0, -1)^T$.

The corresponding NM is shown in Figs. 4(a)–(d). As was to be expected, the damping function $\xi$ was found to differ significantly from zero, but kept a constant (negative) sign starting at $\xi = -0.649$ at $v = 0$ and increasing slowly with $v$. The corresponding invariant manifold, including the modal motion with initial
<figure><img src="image_5.png" /><figcaption>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 $\phi$, (b) $X_2$ versus $v$ and $\phi$, (c) $\Omega$ versus $v$ and $\phi$ and (d) $\xi$ versus $v$ and $\phi$.</figcaption></figure>

conditions $v(0) = v_{\max}$ and $\phi(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 $\rho(T_{\max}) = 0.00962$ with $N_v = 4$ and $N_\phi = 2$.

### 5.1.4. Dissipative system with non-proportional damping

Parameter values: $c_{12} = c_{21} = 0.1$, $c_{ii} = 2\tau\omega_i$ with $\tau = 0.2$. Linear mode under consideration: $\eta = -0.649$, $\omega = 3.177$, $\psi^c = (-0.0324, 0.0001)^T$ and $\psi^s = (0.0043, 0.9994)^T$.

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 $\phi(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 $\rho(T_{\max}) = 0.0103$ with $N_v = 4$ and $N_\phi = 2$.

The invariant manifolds were plotted in Figs. 6(a)–(d) in the phase subspace $(Q_1, Q_2, \dot{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.
<figure><img src="image_6.png" /><figcaption>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 $\phi$, (b) $X_2$ versus v and $\phi$, (c) $\Omega$ versus v and $\phi$ and (d) $\xi$ versus v and $\phi$.</figcaption></figure>

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 $\phi(0) = \phi \in [0, 2\pi]$.

## 5.2. A Van der Pol system

Let us consider the two coupled van der Pol equations

$$ \begin{aligned} \ddot{q}_1 + c_{11}\dot{q}_1 + c_{12}\dot{q}_2 - \varepsilon_1(1 - q_1^2 - \delta_1q_2^2)\dot{q}_1 + q_1 &= 0, \\ \ddot{q}_2 + c_{21}\dot{q}_1 + c_{22}\dot{q}_2 - \varepsilon_2(1 - \delta_2q_1^2 - q_2^2)\dot{q}_2 + 4q_2 &= 0 \end{aligned} \quad (68) $$

with parameter values $\varepsilon_1 = \varepsilon_2 = 0.5$, $\delta_1 = \delta_2 = 1$, $c_{11} = c_{22} = 0$ and $c_{12} = -c_{21} = 0.8$.

Nonlinear motions were computed with $v_{\max} = 2$, $N_v = 4$ and $N_{\phi} = 2$, starting with the modes of the linearized system:

* Mode 1: $\eta = 0.21$, $\omega = 0.89$, $\Psi^c = (-0.011, 0.220)^T$, $\Psi^s = (-0.974, -0.049)^T$,

* Mode 2: $\eta = 0.29$, $\omega = 2.17$, $\Psi^c = (0.046, 0.907)^T$, $\Psi^s = (0.417, -0.021)^T$.
<figure><img src="image_7.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d)</figcaption></figure>

<caption>Fig. 6. Invariant manifold in the phase subspace ($Q_1, Q_2, \dot{Q}_2$): (a) non-gyroscopic conservative system, (b) gyroscopic conservative system, (c) proportional damping and (d) non-proportional damping.</caption>

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.
<figure><img src="image_11.png" /><figcaption>Fig. 7. Frequency ((a), (b)) and damping ((c) and (d)) functions associated with the van der Pol system (68): (a) and (c) nonlinear mode 1 and (b) and (d) nonlinear mode 2.</figcaption></figure>

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 exists 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. [11])

$$ \frac{dv}{d\phi} = v\langle\tau(v)\rangle, \quad \text{where} \quad \langle\tau(\cdot)\rangle = \frac{1}{2\pi} \int_{0}^{2\pi} \frac{\xi(\cdot, \phi)}{\Omega(\cdot, \phi)} d\phi. \qquad (69) $$

Mean damping rates $\langle\tau\rangle$ 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 asymptotically stable ($d/dv\langle\tau\rangle(v_i^*)<0$) 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^*(\phi) = 1.5957 - 0.0059 \cos 2\phi - 0.0010 \cos 4\phi - 0.0430 \sin 2\phi + 0.0004 \sin 4\phi, $$

$$ v_2^*(\phi) = 1.8657 + 0.0461 \cos 2\phi - 0.01761 \cos 4\phi + 0.2146 \sin \phi - 0.0043 \sin 4\phi. $$

giving the limit cycle approximations ($i = 1, 2$)

$$ \left( \begin{array}{c} \mathbf{Q}(t) \\ \dot{\mathbf{Q}}(t) \end{array} \right)_i = v_i^*(\phi(t)) \begin{array}{c} \mathbf{X}_i(v_i^*(\phi(t)), \phi(t)) \\ \mathbf{Y}_i(v_i^*(\phi(t)), \phi(t)) \end{array} \quad \text{with} \quad \dot{\phi} = \Omega_i(v_i^*(\phi), \phi), \quad \phi(0) = \varphi. \qquad (70) $$
<figure><img src="image_12.png" /><figcaption>Fig. 8. Evolution of (τ) versus v for the nonlinear modes 1 (continuous line) and 2 (dashed line).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>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).</figcaption></figure>

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

$$ \frac{d}{dt} \mathbf{Z}(t) = \mathbf{G}(\mathbf{Z}(t)) \quad \text{with} \quad \mathbf{Z} = (\mathbf{Q}^\mathrm{T}, \dot{\mathbf{Q}}^\mathrm{T})^\mathrm{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 $2\pi$-periodic variational linear differential system (see Ref. [11, p. 119])

$$ \frac{d}{d\phi} \Delta \mathbf{Z}(\phi) = \frac{1}{\Omega_i(v_i^*(\phi), \phi)} [\partial_t \mathbf{G}_Z(\mathbf{Z}_i(\phi))] \Delta \mathbf{Z}(\phi). \quad (71) $$

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^*$.

## 6. 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.

## References

[1] R.M. Rosenberg, On nonlinear vibrations of systems with many degree of freedom, *Advances in Applied Mechanics* 242 (9) (1966) 155-242.

[2] A.H. Nayfeh, D.T. Mook, *Nonlinear Oscillations*, Wiley, New York, 1984.

[3] W. Szemplinska-Stupnicka, *The Behavior of Nonlinear Vibrating Systems*, Vols. I and II, Kluwer Academic Publishers, Dordrecht, 1990.

[4] L. Jezequel, C.H. Lamarque, Analysis of nonlinear dynamical systems by normal form theory, *Journal of Sound and Vibration* 149 (3) (1991) 429-459.

[5] S.W. Shaw, C. Pierre, Non-linear normal modes and invariant manifolds, *Journal of Sound and Vibration* 150 (1) (1991) 170-173.

[6] A.F. Vakakis, Nonlinear normal modes (NNMs) and their applications in vibration theory, *Mechanical Systems and Signal Processing* 11 (1) (1997) 3-22.

[7] S.W. Shaw, C. Pierre, Normal modes for nonlinear vibratory systems, *Journal of Sound and Vibration* 164 (1) (1993) 85-124.

[8] E. Pescheck, S.W. Shaw, C. Pierre, A new Galerkin-based approach for accurate nonlinear normal modes through invariant manifolds, *Journal of Sound and Vibration* 249 (5) (2002) 971-993.

[9] S. Bellizzi, R. Bouc, A new formulation for the existence and calculation of non-linear normal modes, *Journal of Sound and Vibration* 287 (2005) 545-569.

[11] J. Hale, *Ordinary Differential Equations*, Wiley-Interscience, New York, 1969.

[12] R. Arquier, F. Pérignon, B. Cochelin, Numerical continuation of nonlinear modes of elastic structures, in: S. Bellizzi, B. Cochelin, C.H. Lamarque (Eds.), *EUROMech 457 Nonlinear Modes of Vibrating Systems*, Press of Ecole Nationale des Travaux Publics de l'Etat, 2004, pp. 17-20, ISBN 2-86 834-119-5.

[13] Z. Xu, Y.K. Chung, Averaging method using generalized harmonic functions for strongly non-linear oscillators, *Journal of Sound and Vibration* 174 (1994) 563-576.