# A new formulation for the existence and calculation of nonlinear normal modes

Sergio Bellizzi, Robert Bouc

► To cite this version:

Sergio Bellizzi, Robert Bouc. A new formulation for the existence and calculation of nonlinear normal modes. Journal of Sound and Vibration, 2005, 287, pp.545-569. 10.1016/j.jsv.2004.11.014. hal-00087990

HAL Id: hal-00087990

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

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
# A new formulation for the existence and calculation of nonlinear normal modes

Sergio Bellizzi*, Robert Bouc

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

A new formulation is presented here for the existence and calculation of nonlinear normal modes in undamped nonlinear autonomous mechanical systems. As in the linear case an expression is developed for the mode in terms of the amplitude, mode shape and frequency, with the distinctive feature that the last two quantities are amplitude and total phase dependent. The dynamic of the periodic response is defined by a one-dimensional nonlinear differential equation governing the total phase motion. The period of the oscillations, depending only on the amplitude, is easily deduced. It is established that the frequency and the mode shape provide the solution to a $2\pi$-periodic nonlinear eigenvalue problem, from which a numerical Galerkin procedure is developed for approximating the nonlinear modes. The procedure is applied to various mechanical systems with two degrees of freedom.

## 1. Introduction

Extending the concept of normal modes of vibrating systems to the case where the restoring forces contain nonlinear terms has been a challenge to many authors, mainly because the principle of linear superposition does not hold for nonlinear systems. However, the existence of “synchronous” periodic oscillations has great potential for applications to nonlinear free and forced vibration problems. The recent review paper by Vakakis [1], discusses the need for the

<footnote>*Corresponding author.</footnote>

<footnote>E-mail address: bellizzi@lma.cnrs-mrs.fr (S. Bellizzi).</footnote>
normal modes approach to be extended to nonlinear theory, and it is concluded that nonlinear normal modes (NNMs) may provide a valuable theoretical tool for understanding some peculiarities of nonlinear systems such as mode bifurcations and nonlinear mode localization.

Following the pioneer work by Rosenberg [2] on conservative systems, several attempts have been made to develop methods of calculating NNMs. These include the harmonic method developed by Szemplinska-Stupnicka [3], the normal form theory [4,5], the invariant manifold method [6,7] (which led to a new definition of NNMs extending the concept to non-conservative systems), the perturbation method [8], the balance harmonic procedure [9], the method of multiple scales [10], and various combinations.

In the present study, a new formulation is presented for the existence and calculation of the synchronous periodic oscillation of an undamped autonomous nonlinear mechanical system. The vector restoring force does not necessarily derive from a potential function. We define motion in line with Rosenberg's definition [2], i.e. a motion where all the material points in the system reaches their extreme values and/or pass through zero simultaneously (NNM). The modal line in the configuration space can be either straight or curved. The NNM is expressed here like a linear normal mode (LNM), in terms of the "amplitude", "phase", "modal vector", and "frequency". The most distinctive feature of this formulation is that the modal vector and frequency are viewed as *amplitude* and *total phase* dependent. The dynamic is defined by a one-dimensional differential equation, governing the total phase motion, from which the period of the oscillations is deduced. The period depends only on the amplitude. The modal vector and the frequency provide the solution to a nonlinear eigenvector-eigenvalue 2π-periodic problem, which makes it possible to calculate these quantities using the classical Galerkin procedure in the space of 2π-periodic functions [10]. It should be noted that the zero-order solution reduces to the Szemplinska-Stupnicka approach [3]. The formulation describes the NNM in terms of synchronous periodic oscillations as well as in terms of a two-dimensional invariant set of the dynamical system.

The case of an odd nonlinear vector restoring force is first considered. A necessary and sufficient condition for the existence of similar NNMs, i.e. modes for which the modal line is a straight line in the configuration space, is given. A basic result as regards the existence and uniqueness of an NNM located outside bifurcating points is also given. The Galerkin computational procedure is applied to two-dimensional mechanical systems in the odd case as well as in the general case. It seems that the presence of internal resonances in the underlying linear system does not require particular attention in these calculations.

## 2. Problem under study

In this study, we present a new formulation for the existence and calculation of synchronous periodic solutions of the undamped autonomous nonlinear n-degrees of freedom (dof) mechanical system

$$[M]\ddot{X}(t) + F(X(t)) = 0, \quad (1)$$

where the mass matrix, $[M]$, is symmetrical and positive definite, and the restoring force vector function, $F(X)$ (including linear and nonlinear terms) is continuously differentiable and such that $F(0) = 0$. The over dots stand for temporal derivatives.
The following assumption will be made throughout this study:

H1. F is a monotonically increasing function, viz.

$$ \forall \mathbf{X}, \forall \mathbf{Y} \in \mathbb{R}^n, \quad (\mathbf{X} - \mathbf{Y})^\mathrm{T}(\mathbf{F}(\mathbf{X}) - \mathbf{F}(\mathbf{Y})) \geq 0. \qquad (2) $$

The linear equation

$$ [\mathbf{M}]\ddot{\mathbf{X}}(t) + [\partial\mathbf{F}_X(\mathbf{0})]\mathbf{X}(t) = \mathbf{0}, \qquad (3) $$

where $[\partial\mathbf{F}_X(\mathbf{0})]$ denotes the Jacobian matrix of F with respect to X, will be called the underlying linear system (or linearized system) associated with the nonlinear system (1). It can be readily verified that: if $\mathbf{X}^*(t)$ is a periodic solution with period $T > 0$ (i.e. $\mathbf{X}^*(t) = \mathbf{X}^*(t+T), \forall t$) of Eq. (1), then $\mathbf{X}^{**}(t) = \mathbf{X}^*(-t)$ is also a T-periodic solution for Eq. (1). Our study will be devoted to the T-periodic solution of Eq. (1) with symmetry $\mathbf{X}^*(t) = \mathbf{X}^*(-t)$. That is the associated Fourier series will contain only cosine terms.

An example of restoring forces associated with geometrical nonlinear behaviour of vibrating beam, plate or shell [11,12] is given by

$$ \mathbf{F}(\mathbf{X}) = \mathbf{L}(\mathbf{X}) + \mathbf{Q}(\mathbf{X}, \mathbf{X}) + \mathbf{C}(\mathbf{X}, \mathbf{X}, \mathbf{X}), \qquad (4) $$

where L, Q and C are linear, bi-linear and tri-linear vector functions. In order that assumption H1 be satisfied the odd part of $\mathbf{F}(\mathbf{X})$ (linear and cubic terms) are of prime importance. In the following, the odd restoring force case will be first considered. The general case will be discussed in Section 4.

### 3. Odd restoring force case

In addition to H1, we assume throughout this section that

$$ \text{H2. } \mathbf{F}(\mathbf{X}) = -\mathbf{F}(-\mathbf{X}). $$

#### 3.1. LNMs

We will first deal briefly with the classical linear case in which $\mathbf{F}(\mathbf{X}) = [\mathbf{K}]\mathbf{X}$, where $[\mathbf{K}]$ is a symmetrical positive definite matrix. The general solution of Eq. (1) can be expressed as

$$ \mathbf{X}(t) = \sum_{l=1}^{n} \psi_l a_l \cos(\omega_l t + \varphi_l), $$

where $(\psi_l, \omega_l^2)$ with $\omega_l > 0$ for $l = 1, \dots, n$ denotes the set of eigenvector-eigenvalue such that

$$ [\mathbf{K}]\psi_l = [\mathbf{M}]\psi_l\omega_l^2, \quad \psi_l^\mathrm{T}[\mathbf{M}]\psi_l = 1 \quad \text{and} \quad \psi_l^\mathrm{T}[\mathbf{M}]\psi_k = 0, \quad l \neq k. $$

The pairs $(a_l, \varphi_l) \in \mathbb{R}_+ \times [0, 2\pi]$ of constant real numbers set the initial conditions. It follows that Eq. (1) possesses n-periodic solutions with period $T_l = 2\pi/\omega_l$, $l = 1, \dots, n$, respectively, which
can be written as (index *l* is omitted)

$$X(t) = \psi a \cos \Phi(t), \qquad (5)$$

$$\dot{\Phi}(t) = \omega, \quad \Phi(0) = \varphi. \qquad (6)$$

The scalar function $\Phi(\cdot)$ describes the dynamic operating on the modal line (a straight line) in the configuration space.

In cases where the matrix $[K]$ is not symmetrical but only positive the same result holds. The eigenvectors will no longer be $[M]$-orthogonal, however, in this case.

## 3.2. NNMs

As suggested by the linear case, a synchronous periodic solution to Eq. (1) is sought in the form

$$X(t) = \Psi(a, \Phi(t))a \cos \Phi(t), \qquad (7)$$

$$\dot{\Phi}(t) = \Omega(a, \Phi(t)), \quad \Phi(0) = \varphi. \qquad (8)$$

Eq. (7) describes the geometrical behaviour of the solution in the configuration space and the dynamics of the periodic response is defined by the one-dimensional differential equation, Eq. (8), governing the total phase motion $\Phi(\cdot)$.

Four quantities are involved in Eqs. (7) and (8): two given scalar quantities, the *amplitude* ($a \ge 0$) and the *phase* ($\varphi \in [0, 2\pi]$) which set the initial conditions of the vibration in NNM motion on the configuration space. Two unknown functions, $\Psi$ and $\Omega$, that will be referred to here as the *modal vector* and the *modal or resonance frequency* of the NNM, respectively. To ensure that parameter $a$ appropriately characterizes the amplitude of the modal line in Eq. (7), a normalization condition on the modal vector $\Psi$ is required. In this study, we will adopt (without lost of generality) the following condition:

$$\Psi^T(a, \Phi)[M]\Psi(a, \Phi) = 1. \qquad (9)$$

The modal vector $\Psi$ and the frequency $\Omega$ will be viewed here as *amplitude* and *total phase dependent*, and denoted $\Psi = \Psi(a, \Phi)$, $\Omega = \Omega(a, \Phi)$, and will be searched for as $2\pi$-periodic functions with respect to $\Phi$ for fixed $a$. More specifically, according to the symmetry condition H2, these functions will be searched for as *even periodic functions* with respect to $\Phi$ with period $2\pi$, that is with Fourier series

$$\Omega(a, \Phi) = \sum_{k=0}^{\infty} \Omega_{2k}(a) \cos 2k\Phi, \quad \Psi(a, \Phi) = \sum_{k=0}^{\infty} \Psi_{2k}(a) \cos 2k\Phi.$$

It follows that the Fourier series of $X(t)$ in terms of $\Phi(t)$ will contain only odd cosine terms. From Eq. (8), it results that

$$dt = \frac{1}{\Omega(a, \Phi)} d\Phi.$$
Once the (positive) frequency has been determined, the period of the oscillations will be given by

$$T(a) = \int_{\varphi}^{\varphi+2\pi} \frac{1}{\Omega(a, \varphi)} d\varphi = 2 \int_{0}^{\pi} \frac{1}{\Omega(a, \varphi)} d\varphi, \quad (10)$$

and the result depends only on $a$. It can be readily verified that $\Phi(t+T(a)) = \Phi(t)+2\pi$ and $\Phi(t+T(a)/2) = \Phi(t)+\pi$ yielding $\mathbf{X}(t) = \mathbf{X}(t+T(a))$ and $\mathbf{X}(t) = -\mathbf{X}(t+T(a)/2)$.

Finally, Eq. (7) together with

$$\dot{\mathbf{X}}(t) = a\Omega(a, \Phi(t)) \left( \frac{\partial \mathbf{\Psi}(a, \Phi(t))}{\partial \Phi} \cos \Phi(t) - \mathbf{\Psi}(a, \Phi(t)) \sin \Phi(t) \right), \quad (11)$$

define a “synchronous” periodic oscillation [2], where all the material points in the system pass through zero ($\mathbf{X}(t^*) = \mathbf{0}$, $\Phi(t^*) = \pi/2$) or reach their extreme values ($\mathbf{X}(t^{**}) = \mathbf{0}$, $\Phi(t^{**}) = 2\pi$) simultaneously. The set, in the phase space, defined by Eqs. (7) and (11) taking the initial conditions $a$ and $\varphi$ to be independent variables is an invariant set for the dynamical system associated with the equation of motion, Eq. (1). Thus, this formulation gives also a characterization of the NNM in the framework of invariant manifold [6,7].

The proposed formulation differs from that given in Ref. [3]

$$\mathbf{X}(t) = a\mathbf{\Psi}(a) \cos(\Omega(a)t + \varphi), \quad (12)$$

where the frequency and the mode shape are only dependent on the amplitude. Taking into account the total phase in the definition of the modal vector allows us considering the modal lines (in the configuration space) either straight or curved.

### 3.3. A periodic eigenvalue–eigenfunction problem

The objective now is to characterize the pair $(\Omega, \mathbf{\Psi})$ which defines an NNM. Substituting Eq. (7) into Eq. (1), using Eq. (8), we get

$$[\mathbf{M}] \frac{\mathrm{d}\left(\frac{\mathrm{d}(\mathbf{\Psi}(a, \Phi) \cos \Phi)}{\mathrm{d}\Phi} \Omega(a, \Phi)\right)}{\mathrm{d}\Phi} \Omega(a, \Phi) + \frac{1}{a} \mathbf{F}(a\mathbf{\Psi}(a, \Phi) \cos \Phi) = \mathbf{0}, \quad (13)$$

where $\Phi$ looks like an independent variable. This variable will therefore be denoted $\phi$ in what follows.

Differentiating with respect to $\phi$, using the differential rule $2\Omega\partial\Omega/\partial\phi = \partial\Omega^2/\partial\phi$, Eq. (13) becomes

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

where the differential operator $\mathbf{L}$ is given by

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

The frequency function, $\Omega$, appears in Eq. (14) by its square, $\Omega^2$. So, in the following, we will look for the pair $(\Omega^2, \mathbf{\Psi})$.
Recalling

$$ \Psi^T[M]\Psi = 1, \qquad (16) $$

we must solve, for the fixed parameter $a$, the $2\pi$-periodic nonlinear eigenvalue $(\Omega^2)$-eigenvector $(\Psi)$ problem defined by Eqs. (14) and (16). These equations reduce to $[K]\Psi = [M]\Psi\Omega^2$ in the linear case.

We will focus on periodic solutions $(\Omega^2(a, .), \Psi(a, .))$ which satisfy the following properties:

* $\phi \to (\Omega^2(a, \phi), \Psi(a, \phi))$ are real-valued functions and (one time, two times) continuously differentiable on $[0, 2\pi]$,

* $\forall \phi, (\Omega^2(a, \phi), \Psi(a, \phi)) = (\Omega^2(a, -\phi), \Psi(a, -\phi)) = (\Omega^2(a, \phi + \pi), \Psi(a, \phi + \pi))$.

Such solutions will be named *well-defined* periodic solutions.

It can be shown (see Appendix A) that for a well-defined periodic solution of Eqs. (14) and (16) we always have

$$ \forall (a, \phi), \Omega^2(a, \phi) > 0, \qquad (17) $$

and $\Omega^2$ can be expressed in terms of $\Psi$ as

$$ \Omega^2(a, \phi) = \frac{2}{a \sin^2 \phi} \left( \int_0^\phi \exp\left(-\int_\sigma^\phi \beta(a, v) \, dv\right) I(a, \sigma) \, d\sigma \right), \qquad (18) $$

where

$$ I(a, \phi) = \sin \phi \Psi^T(a, \phi) F(\Psi(a, \phi)) a \cos \phi, $$

and

$$ \beta(a, \phi) = 2\gamma^2(a, \phi) \frac{\cos \phi}{\sin \phi} \quad \text{with} \quad \gamma^2(a, \phi) = \frac{\partial \Psi^T(a, \phi)}{\partial \phi} [M] \frac{\partial \Psi(a, \phi)}{\partial \phi}. $$

Note that inequality (17) ensures the existence of a real positive resonance frequency function $\Omega$ and Eq. (18) reduces to $\Omega^2 = \Psi^T[K]\Psi$ in the linear case.

Moreover, if we assume that the eigenvalues associated with the pair of matrices $[M], [\partial_X F(0)]$ are distinct and without resonance relations between the associated frequencies then we can prove that for each $a$ in some neighbourhood of $a = 0$, there exist $n$ well-defined solutions to Eqs. (14) and (16). Each solution is unique in some neighbourhood of $a = 0$ and $(\omega_l^2, \psi_l)$ where $(\omega_l^2, \psi_l)$ denotes the corresponding normal mode of the underlying linear system (3).

Furthermore, as $a \to 0$, the following limits hold:

$$ \sup_{\phi \in [0, 2\pi]} |\Omega^2(a, \phi) - \omega_l^2| \to 0, \quad \sup_{\phi \in [0, 2\pi]} \|\Psi(a, \phi) - \psi_l\| \to 0, $$

$$ \sup_{\phi \in [0, 2\pi]} \left| \sin \phi \frac{\partial \Omega^2}{\partial \phi} \right| \to 0, \quad \sup_{\phi \in [0, 2\pi]} \left| \frac{\partial \Psi}{\partial \phi} \right| \to 0, \quad \sup_{\phi \in [0, 2\pi]} \left| \frac{\partial^2 \Psi}{\partial \phi^2} \right| \to 0. $$

Consequently there exist $n$ NNMs which can be viewed as a continuation of the LNM of the underlying linear system (3).
### 3.4. A necessary and sufficient condition for the existence of similar NNMs (SNNMs)

As described in Ref. [2], a nonlinear mode is “similar” if the motion in the configuration space is such that

$$ \frac{X_i(t)}{X_j(t)} = (\text{Constant})_{ij}, \qquad (19) $$

where $X_i$ denotes the $i$th component of $\mathbf{X}$.

This leads to the following statement: an NNM is similar if and only if the $[\mathbf{M}]$-normalized vector $\mathbf{\Psi}$ does not depend on $\mathbf{\Phi}$.

In this case, Eq. (14), using Eq. (16), reduces to

$$ \Omega^2 \cos \phi + \frac{1}{2} \sin \phi \frac{\partial \Omega^2}{\partial \phi} = \mathbf{\Psi}^\mathrm{T} \frac{1}{a} \mathbf{F}(a\mathbf{\Psi} \cos \phi), \qquad (20) $$

and substituting Eq. (20) into Eq. (14) yields

$$ \mathbf{F}(a\mathbf{\Psi} \cos \phi) = (\mathbf{\Psi}^\mathrm{T} \mathbf{F}(a\mathbf{\Psi} \cos \phi))[\mathbf{M}]\mathbf{\Psi}. \qquad (21) $$

This leads to the following statement: the nonlinear system (1) admits an SNNM if and only if Eq. (21) possesses a $[\mathbf{M}]$-normalized $\phi$-independent vector solution $\mathbf{\Psi}$. The associated frequency is given by

$$ \Omega^2(a, \phi) = \frac{2}{a \sin^2 \phi} \int_0^\phi \mathbf{\Psi}^\mathrm{T} \mathbf{F}(a\mathbf{\Psi} \cos \sigma) \sin \sigma \, d\sigma. \qquad (22) $$

This expression defines a regular, even, $2\pi$-periodic function with respect to $\phi$ with a Fourier series of the form $\Omega^2(a, \phi) = \sum_{k=0}^\infty \Omega_{2k}^2(a) \cos 2k\phi$.

*Comment:* It can be easily checked that any function of the following form yields SNNMs:

* $\mathbf{F}(\mathbf{X}) = [\mathbf{K}]\mathbf{X} + f(\mathbf{X})[\mathbf{K}_1]\mathbf{X}$ where $[\mathbf{K}]$ and $[\mathbf{K}_1]$ are symmetrical and definite positive matrices such that the eigenvalue problem associated with the pair of matrices ($[\mathbf{M}]$, $[\mathbf{K}]$) and ($[\mathbf{M}]$, $[\mathbf{K}_1]$) leads to the same eigenvectors, and $f(\mathbf{X})$ is some even, positive, smooth scalar function (in that case the mode shapes of the SNNMs coincide with those of the underlying linear system (see example below);

* $\mathbf{F}(\mathbf{X}) = \alpha[\mathbf{M}]\mathbf{X} + \mathbf{G}(\mathbf{X})$ where $\alpha \ge 0$ and $G_i(\mathbf{X})$, $i = 1, 2, \dots, n$, are homogeneous polynomials of degree $2k + 1$.

**Example (a multidimensional Duffing equation).** Let

$$ \mathbf{F}(\mathbf{X}) = [\mathbf{K}]\mathbf{X} + \lambda(\mathbf{X}^\mathrm{T}[\mathbf{C}]\mathbf{X})[\mathbf{K}]\mathbf{X}, $$

where $[\mathbf{K}]$ and $[\mathbf{C}]$ are real square matrices, $[\mathbf{K}]$ is symmetrical and positive definite and $\lambda$ is a real number. Eq. (21) reduces to

$$ [\mathbf{K}]\mathbf{\Psi} = (\mathbf{\Psi}^\mathrm{T}[\mathbf{K}]\mathbf{\Psi})[\mathbf{M}]\mathbf{\Psi}. $$
Hence the mode shapes of the SNNMs coincide with those of the underlying linear system ($\lambda = 0$). Let $(\omega_p, \psi_p)$ be a normal mode of the linear system. Using $\psi_p^T[M]\psi_p = 1$ and $\psi_p^T[K]\psi_p = \omega_p^2$, Eq. (22) gives

$$ \Omega_p^2(a, \phi) = \omega_p^2 + \lambda \frac{\omega_p^2 a^2}{4} (\psi_p^T[C]\psi_p)(3 + \cos 2\phi). $$

The pair $(\Omega_p^2, \psi_p)$ completely defines an SNNM which tends towards $(\omega_p^2, \psi_p)$ as $a$ tends towards zero. The normal mode can be a softening or hardening mode depending on the sign of $\lambda(\psi_p^T[C]\psi_p)$ (provided hypothesis H1 is fulfilled). The corresponding period of the oscillations is given by (10)

$$ T(a) = \frac{4}{\omega_p \sqrt{1 + \lambda \psi_p^T[C] \psi_p a^2}} K \frac{\lambda (\psi_p^T[C] \psi_p) a^2}{2(1 + \lambda (\psi_p^T[C] \psi_p) a^2)}, $$

where $K(.)$ denotes the elliptic function of the first kind. Note that $T(a = 0) = 2\pi/\omega_p$.

It is worth taking this example to compare our approach with the normal form approach as described, for instance, in Ref. [5]. The two approaches both yield the same geometry of the manifold in the phase space and the same oscillation period. The main difference focuses on the representation of the dynamic. The two approaches are of course equivalent.

### 3.5. A Galerkin procedure for the calculation of the NNMs

Now it is proposed to obtain accurate approximate solutions to Eqs. (14) and (16). We have to solve a periodic differential equation with known period ($2\pi$). To do this, a Galerkin method is implemented [10]. The unknown functions (square-frequency $\Omega^2$ and mode shape $\Psi$), are expanded into a finite Fourier series with respect to the variable $\phi$ according to

$$ \Omega_m^2(a, \phi) = \sum_{k=0}^{m} \Omega_{m,2k}^2(a) \cos 2k\phi, \quad \Psi_m(a, \phi) = \sum_{k=0}^{m} \Psi_{m,2k}(a) \cos 2k\phi. $$

The determining equations read, for $k = 0, \dots, m$

$$ \int_0^{2\pi} \left( \mathbf{L}(\Omega_m^2, \mathbf{\Psi}_m; \phi) + \frac{1}{a} \mathbf{F}(a\mathbf{\Psi}_m \cos \phi) - [\mathbf{M}]\mathbf{\Psi}_m \left( \Omega_m^2 \cos \phi + \frac{1}{2} \frac{\partial \Omega_m^2}{\partial \phi} \sin \phi \right) \right) \times \cos(2k+1)\phi \, d\phi = \mathbf{0}, \quad (23) $$

$$ \int_0^{2\pi} (\mathbf{\Psi}_m^T[M]\mathbf{\Psi}_m - 1) \cos(2k) \phi \, d\phi = 0. \quad (24) $$

This constitutes a set of $(m+1)n + m + 1$ nonlinear equations for the $(m+1)n + m + 1$ unknown coefficients $\Omega_{m,2k}^2$ and $\Psi_{m,2k}$. It should be noted that these equations are in agreement with the set of functions usually used with this method, because due to the symmetrical properties of the
various terms with respect to the variable $\phi$, the following equations

$$ \int_{0}^{2\pi} \left( \mathbf{L}(\Omega_{m}^{2}, \mathbf{\Psi}_{m}; \phi) + \frac{1}{a} \mathbf{F}(a\mathbf{\Psi}_{m} \cos \phi) - [\mathbf{M}]\mathbf{\Psi}_{m} \left( \Omega_{m}^{2} \cos \phi + \frac{1}{2} \frac{\partial \Omega_{m}^{2}}{\partial \phi} \sin \phi \right) \right) \times \cos 2k\phi \, d\phi = \mathbf{0}, $$

$$ \int_{0}^{2\pi} (\mathbf{\Psi}_{m}^{\mathrm{T}}[\mathbf{M}]\mathbf{\Psi}_{m} - 1) \cos(2k + 1)\phi \, d\phi = 0 $$

are always satisfied for $k = 1, \dots, m$.

For given $a$ and $m$, Eqs. (23) and (24) have to be solved numerically. A Newton-Raphson method has been implemented with exact Jacobian matrices, together with an incremental-continuation procedure with respect to the parameters $a$ and/or $m$. For $a = a + \Delta a$ or for $m = m + 1$, the previous approximate values are used as the starting point. When $m$ is increased, the $(m+1)$-order coefficients are initially taken to be equal to zero. Solutions of the reduced equations for $m = 0$,

$$ \int_{0}^{2\pi} \frac{1}{a} \mathbf{F}(a\mathbf{\Psi}_{0,0} \cos \phi) \cos \phi \, d\phi - 2\pi[\mathbf{M}]\mathbf{\Psi}_{0,0}\Omega_{0,0}^{2} = \mathbf{0}, $$

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

can be used to start the procedure, especially with small values of $a$, where the solution is close to the corresponding LNM.

The accuracy of the approximation can be checked, evaluating with respect to $m$, the $L^2$-norm of the residual

$$ r_m(a) = \int_0^{2\pi} \left\| \mathbf{L}(\Omega_m^2, \mathbf{\Psi}_m; \phi) + \frac{1}{a} \mathbf{F}(a\mathbf{\Psi}_m \cos \phi) - [\mathbf{M}]\mathbf{\Psi}_m \left( \Omega_m^2 \cos \phi + \frac{1}{2} \frac{\partial \Omega_m^2}{\partial \phi} \sin \phi \right) \right\|^2 d\phi \\ + \int_0^{2\pi} (\mathbf{\Psi}_m^\mathrm{T}[\mathbf{M}]\mathbf{\Psi}_m - 1)^2 d\phi, \qquad (25) $$

and/or this sup-norm.

Once the $\Omega_{m,2k}^2$'s and $\mathbf{\Psi}_{m,2k}$'s have reached the desired level of accuracy at the desired amplitude level $a$, the differential equation (8) can be solved numerically. The resulting time history, $\Phi(t)$, allows us to derive the corresponding time histories of the displacement (7) and velocity (11).

## 3.6. Examples

### 3.6.1. Example 1

The first example is taken from Ref. [7]

$$ \ddot{x}_1 + \omega_1^2 x_1 + 0.405 x_1^3 + 1.34 x_1^2 x_2 + 1.51 x_1 x_2^2 + 0.349 x_2^3 = 0, \qquad (26) $$

$$ \ddot{x}_2 + \omega_2^2 x_2 + 0.448 x_1^3 + 1.51 x_1^2 x_2 + 1.05 x_1 x_2^2 + 4.580 x_2^3 = 0, \qquad (27) $$

where $\omega_1 = 0.689$ and $\omega_2 = 3.244$ denote the natural frequencies of the underlying linear system.
The first NNM is depicted in Figs. 1–4 as a continuation of the linearized mode $\Omega_{0,0}^2 = 0.689^2$, $\Psi_{0,0} = (1,0)^T$. The residual $r_m(a)$ is shown in the Log-linear plot in Fig. 1 for various values of $m$. The amplitude $a$ is limited to the domain $[0, 1.88]$. The residual decreases with $m$ but the more $a$ increases, the less the residual decreases. At greater amplitudes the Newton–Raphson method fails to converge, which indicates that a bifurcation is present. A similar case will be considered in Example 2. The mode shape $\Psi(a, \phi)$ and the frequency $\Omega(a, \phi)$ are shown in Fig. 2(a)–(c) for a ($m=6$)-harmonics Galerkin approximation. These surfaces are similar to that of the linear normal mode at small values of $a$. When $a$ increases, the fluctuations in $\Psi(a, 0)$ go beyond the range of small motions, and large-amplitude motion is generated. At fixed $\phi$, the frequency $\Omega$ increases with the amplitude $a$. The maximum values of $\Omega$ in the cross-sections defined by $a = C'$ are obtained at $\phi = 0$ and $\pi$ by periodicity. The associated period, which depends only on the amplitude $a$, is shown in (d). The hardening behaviour of the system is confirmed.

Fig. 3 gives a picture of the invariant manifold in the phase subspace ($X_1, X_2, Y_1$), where the axis $Y_1$ corresponds to the velocity variable $\dot{x}_1$. This mode is non-similar and it can be established that this system does not have a similar mode.

The time histories (over two periods) of the periodic motion corresponding to the amplitude level $a = 1.88$ are compared with numerical simulation in Fig. 4. On the one hand, Eq. (8) is numerically solved with the initial condition $\varphi = 0$ and Eq. (7) is performed using the frequency $\Omega(1.88, \Phi(t))$ and mode shape $\Psi(1.88, \Phi(t))$ with ($m=6$)-harmonics, while on the other hand, a direct numerical integration of Eq. (1) with the corresponding initial conditions ($1.88\Psi_1(1.88,0), 0, 1.88\Psi_2(1.88,0), 0$) in the phase space ($X_1, Y_1, X_2, Y_2$) is carried out. In both cases, the Runge–Kutta method was used to solve the associated differential equations. The two solutions are indistinguishable. The behaviour of the total phase is quasi-linear versus time, due to the small fluctuations of the modal frequency with respect to $\phi$ (Fig. 4(a)).

Similar results were obtained for the second NNM as a continuation of the corresponding linearized mode $\Omega_{0,0}^2 = 3.244^2$, $\Psi_{0,0} = (0,1)$. The $m=3$ order approximation was sufficiently accurate to describe its behaviour.

### 3.6.2. Example 2

We now consider the two dof nonlinear system

$$ \ddot{x}_1 + x_1 + (x_1 - x_2)^3 = 0, \quad (28) $$

$$ \ddot{x}_2 - (x_1 - x_2)^3 = 0. \quad (29) $$

We first focus on the NNM, starting from $\Omega_{0,0}^2 = 0$, $\Psi_{0,0} = (0,1)^T$, in the range $a \in [0, 0.5]$. The convergence of the method is illustrated in Fig. 5. The continuation method fails to converge for $a \approx 0.27$ indicating that a bifurcation is present. A more sophisticated continuation method, such as the asymptotic numerical method, could be used to locate the various branches. The ($m=6$)-order Galerkin approximation is sufficiently accurate to describe the NNM, except in a narrowband around $a = 0.27$ which is omitted in Figs. 5–7. However, the calculated branch for $a > 0.27$ seems again to be the continuation of the linear mode.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Example 1—residual (Eq. (25)) versus amplitude *a*: *m* = 0 (continuous line), *m* = 1 (small dashed line), *m* = 5 (large dashed line) and *m* = 6 (small-large dashed line).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Example 1—mode shape, frequency and period. In plots (a)–(c): *m* = 6. In plot (d): *m* = 0 (continuous line), *m* = 1, 5 and 6 (dashed lines).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Example 1—invariant manifold in the phase subspace ($X_1, Y_1, X_2$) with ($m = 6$).</figcaption></figure>

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

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

<figure><img src="image_7.png" /><figcaption>Fig. 4. Example 1—time histories of the periodic motion: $m = 6$, $a = 1.88$ and $\varphi = 0$. In plot (a): $\Phi$ versus time obtained by solving numerically Eq. (8). In plots (b) and (c): $x_1$ and $x_2$ versus time obtained from Eqs. (7) and (8) (dashed line) and by direct simulation of Eqs. (6) and (27) (continuous line).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 5. Example 2—residual (Eq. (25)) versus amplitude *a*: *m* = 0 (continuous line), *m* = 1 (small dashed line), *m* = 5 (large dashed line) and *m* = 6 (small-large dashed line).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Example 2—mode shape, frequency and period. In plots (a)–(c): *m* = 6. In plot (d): *m* = 0 (continuous line), *m* = 1, 5 and 6 (dashed line).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 7. Example 2—invariant manifold in the phase space ($X_2, Y_2, X_1$) with ($m=6$).</figcaption></figure>

The following ($m=6$)-order approximations were obtained for $a = 0.5$:

$$ \Omega^2(0.5, \phi) = 0.10894 + 0.01249 \cos 2\phi + 0.00422 \cos 4\phi + 0.00371 \cos 6\phi \\ - 0.0000 \cos 8\phi - 0.00019 \cos 10\phi - 0.00006 \cos 12\phi, $$

$$ \Psi_1(0.5, \phi) = 0.00746 + 0.23665 \cos 2\phi + 0.01434 \cos 4\phi - 0.00240 \cos 6\phi \\ - 0.00120 \cos 8\phi - 0.00005 \cos 10\phi + 0.00008 \cos 12\phi, $$

$$ \Psi_2(0.5, \phi) = 0.98576 - 0.00353 \cos 2\phi - 0.01402 \cos 4\phi - 0.00158 \cos 6\phi \\ + 0.00019 \cos 8\phi + 0.00014 \cos 10\phi + 0.00001 \cos 12\phi. $$

It is worth noting that the series expansions are slowly decreasing.

The following truncated expansions have been obtained for the second NNM:

$$ \Omega^2(0.5, \phi) = 1.35188 + 0.13969 \cos 2\phi + 0.00160 \cos 4\phi - 0.00002 \cos 6\phi \\ + 0.00000 \cos 8\phi - 0.00000 \cos 10\phi - 0.0000 \cos 12\phi, $$

$$ \Psi_1(0.5, \phi) = 0.97006 - 0.00340 \cos 2\phi + 0.00002 \cos 4\phi + 0.00000 \cos 6\phi \\ - 0.00000 \cos 8\phi + 0.00000 \cos 10\phi - 0.00000 \cos 12\phi, $$

$$ \Psi_2(0.5, \phi) = -0.24263 - 0.0133 \cos 2\phi + 0.00029 \cos 4\phi - 0.00000 \cos 6\phi \\ + 0.00000 \cos 8\phi - 0.00000 \cos 10\phi - 0.00000 \cos 12\phi $$

as a continuation of the LNM $\Omega^2_{0,0} = 1$, $\Psi_{0,0} = (1,0)^T$.

The two modes are obviously not affected in the same way by the nonlinear terms. The similarity property is almost satisfied by the second mode whereas the first nonlinear mode is strongly non-similar.
## 4. General restoring force case

Throughout this section, we only assume H1.

### 4.1. NNMs

In the case where $F(X)$ is not necessarily an odd function, the NNMs take the following form:

$$X(t) = \Psi(a, \Phi(t))a \cos \Phi(t) + a\mathbf{B}(a, \Phi(t)), \quad (30)$$

$$\dot{\Phi}(t) = \Omega(a, \Phi(t)), \quad \Phi(0) = \varphi, \quad (31)$$

where $\mathbf{B}$ denotes a vector function.

Five quantities are involved in Eqs. (30) and (31): two given scalar quantities, the amplitude ($a>0$) and the phase ($\varphi \in [0, 2\pi]$), which set the initial conditions of the vibration in NNM motion on the configuration space. Three unknown functions, $\Psi$, $\Omega$ and $\mathbf{B}$ that will be referred to here as the modal vector, the modal or resonance frequency and the bias term of the NNM, respectively. The bias term has been added to balance the even cosine terms if the nonlinear function includes even terms. These functions will be searched for as even periodic function with respect to $\Phi$ with period $\pi$, that is with Fourier series

$$\mathbf{B}(a, \varphi) = \sum_{k=0}^{\infty} \mathbf{B}_{2k}(a) \cos 2k\varphi.$$

As in the odd restoring force case, the proposed formulation (30) differs from that given in Ref. [3] where an approximate solution of the form

$$X(t) = a\Psi(a) \cos(\Omega(a)t + \varphi) + a\mathbf{B}(a) \quad (32)$$

has been considered, ignoring the phase dependence on the frequency, mode shape and bias term.

### 4.2. A periodic eigenvalue–eigenfunction problem

The objective now is to characterize the three functions $(\Omega, \Psi, \mathbf{B})$ which define an NNM. Substituting Eq. (30) into Eq. (1), Eqs. (14) and (16) become

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

$$\Psi^T[\mathbf{M}]\Psi = 1, \quad (34)$$

where the operator $\mathbf{L}$ is still given by Eq. (15).

Apparently we have to solve two equations for three unknowns. Nevertheless we look for symmetrical $\pi$-periodic solutions $(\Omega^2(a, .), \Psi(a..), \mathbf{B}(a, .))$ which satisfy the following properties
(well-defined periodic solutions):

* $\phi \to (\Omega^2(a, \phi), \Psi(a, \phi), \mathbf{B}(a, \phi))$ are real-valued functions and (one time, two times)-continuously differentiable on $[0, 2\pi]$,

* $\Omega^2(a, \phi) = \Omega^2(a, -\phi) = \Omega^2(a, \phi + \pi)$, $\Psi(a, \phi) = \Psi(a, -\phi) = \Psi(a, \phi + \pi)$, and $\mathbf{B}(a, \phi) = \mathbf{B}(0, -\phi) = \mathbf{B}(a, \phi + \pi)$.

Consequently, using the decomposition

$$F(\Psi(a, \phi)a \cos \phi + a\mathbf{B}(a, \phi)) = F_{\text{odd}}(a, \Psi, \mathbf{B}; \phi) + F_{\text{even}}(a, \Psi, \mathbf{B}; \phi),$$

where $F_{\text{odd}}(a, \Psi, \mathbf{B}; \phi)$ (resp. $F_{\text{even}}(a, \Psi, \mathbf{B}; \phi)$) stands for the odd (resp. even) cosine terms in the Fourier series of $F(\Psi(a, \phi)a \cos \phi + a\mathbf{B}(a, \phi))$ with respect to $\phi$, the above problem (defined by Eqs. (33) and (34)) can be rewritten into

$$L(\Omega^2, \Psi; \phi) + \frac{1}{a} F_{\text{odd}}(a, \Psi, \mathbf{B}; \phi) - [M]\Psi(\Omega^2 \cos \phi + \frac{1}{2} \frac{\partial \Omega^2}{\partial \phi} \sin \phi) = 0, \quad (35)$$

$$\Omega^2[\mathbf{M}] \frac{\partial^2 \mathbf{B}}{\partial \phi^2} + \frac{1}{2} \frac{\partial \Omega^2}{\partial \phi} [\mathbf{M}] \frac{\partial \mathbf{B}}{\partial \phi} + \frac{1}{a} F_{\text{even}}(a, \Psi, \mathbf{B}; \phi) = 0, \quad (36)$$

$$\Psi^T[\mathbf{M}]\Psi = 1. \quad (37)$$

In the case of a restoring vector force function of the form (4)

$$F(\mathbf{X}) = [\mathbf{K}]\mathbf{X} + \mathbf{Q}(\mathbf{X}, \mathbf{X}) + \mathbf{C}(\mathbf{X}, \mathbf{X}, \mathbf{X}), \quad (38)$$

it is easy to show that

$$\begin{aligned} F_{\text{even}}(a, \Psi, \mathbf{B}; \phi) = [\mathbf{K}]a\mathbf{B} &+ \mathbf{Q}(a\Psi \cos \phi, a\Psi \cos \phi) \\ &+ \mathbf{Q}(a\mathbf{B}, a\mathbf{B}) + \mathbf{C}(a\Psi \cos \phi, a\Psi \cos \phi, a\mathbf{B}) \\ &+ \mathbf{C}(a\Psi \cos \phi, a\mathbf{B}, a\Psi \cos \phi) \\ &+ \mathbf{C}(a\mathbf{B}, a\Psi \cos \phi, a\Psi \cos \phi) + \mathbf{C}(a\mathbf{B}, a\mathbf{B}, a\mathbf{B}), \end{aligned}$$

and

$$\begin{aligned} F_{\text{odd}}(a, \Psi, \mathbf{B}; \phi) = [\mathbf{K}]a\Psi \cos \phi &+ \mathbf{Q}(a\Psi \cos \phi, a\mathbf{B}) \\ &+ \mathbf{Q}(a\mathbf{B}, a\Psi \cos \phi) + \mathbf{C}(a\Psi \cos \phi, a\Psi \cos \phi, a\Psi \cos \phi) \\ &+ \mathbf{C}(a\Psi \cos \phi, a\mathbf{B}, a\mathbf{B}) + \mathbf{C}(a\mathbf{B}, a\Psi \cos \phi, a\mathbf{B}) \\ &+ \mathbf{C}(a\mathbf{B}, a\mathbf{B}, a\Psi \cos \phi), \end{aligned}$$

where the last equation can be re-written as

$$F_{\text{odd}}(a, \Psi, \mathbf{B}; \phi) = a \cos \phi [\partial_x F(a\mathbf{B})] \Psi + a^3 \cos^3 \phi C(\Psi, \Psi, \Psi). \quad (39)$$

Substituting Eq. (39) into Eq. (35) and using the same arguments as in Appendix A, it can be shown that a well-defined periodic solution of Eqs. (35)–(37) satisfies

$$\forall (a, \phi), \quad \Omega^2(a, \phi) > 0. \quad (40)$$
### 4.3. A Galerkin procedure for the calculation of the NNMs

As in the odd restoring force case, an accurate approximate solutions of Eqs. (35)–(37) can be obtained implementing a Galerkin method. The unknown functions are searched for in the form

$$ \Omega_m^2(a, \phi) = \sum_{k=0}^{m} \Omega_{m,2k}^2(a) \cos 2k\phi, \quad \Psi_m(a, \phi) = \sum_{k=0}^{m} \Psi_{m,2k}(a) \cos 2k\phi, $$

$$ \mathbf{B}_m(a, \phi) = \sum_{k=0}^{m} \mathbf{B}_{m,2k}(a) \cos 2k\phi, $$

where *m* denotes the order of the truncated series.

The corresponding determining equations read

$$ \int_{0}^{2\pi} \left( \mathbf{L}(\Omega_{m}^{2}, \mathbf{\Psi}_{m}; \phi) + \frac{1}{a} \mathbf{F}_{\text{odd}}(a, \mathbf{\Psi}_{m}, \mathbf{B}_{m}; \phi) - [\mathbf{M}]\mathbf{\Psi}_{m} \left( \Omega_{m}^{2} \cos \phi + \frac{1}{2} \frac{\partial \Omega_{m}^{2}}{\partial \phi} \sin \phi \right) \right) \times \cos(2k+1)\phi \, d\phi = \mathbf{0}, $$

$$ \int_{0}^{2\pi} \left( \Omega_{m}^{2} [\mathbf{M}] \frac{\partial^{2}\mathbf{B}_{m}}{\partial\phi^{2}} + \frac{1}{2} \frac{\partial\Omega_{m}^{2}}{\partial\phi} [\mathbf{M}] \frac{\partial\mathbf{B}_{m}}{\partial\phi} + \frac{1}{a} \mathbf{F}_{\text{even}}(a, \mathbf{\Psi}_{m}, \mathbf{B}_{m}; \phi) \right) \cos 2k\phi \, d\phi = \mathbf{0}, $$

$$ \int_{0}^{2\pi} (\mathbf{\Psi}_{m}^{\mathrm{T}}[\mathbf{M}]\mathbf{\Psi}_{m} - 1) \cos 2k\phi \, d\phi = 0, $$

for $k = 0, \dots, m$. This constitutes a set of $2n(m+1)+m+1$ nonlinear equations for the $2n(m+1)+m+1$ unknown coefficients $\mathbf{\Psi}_{m,2k}, \mathbf{B}_{m,2k}$, and $\Omega_{m,2k}^2$.

As previously, an incremental-continuation numerical procedure based on the Newton-Raphson algorithm can be implemented.

### 4.4. Example 3

Here we consider a two dof nonlinear system (see Ref. [13]) composed of a mass *m* connected to two springs (see Fig. 8). Under the assumption of large displacement the strain energy is given by

$$ W = \frac{1}{2} k_1 e_1^2 + \frac{1}{2} k_2 e_2^2, $$

where

$$ e_i = \frac{x_i}{L} + \frac{1}{2} \left( \left( \frac{x_1}{L} \right)^2 + \left( \frac{x_2}{L} \right)^2 \right), $$

denotes the strain-displacement of the *r*th spring, the equations of motion can be written, assuming $L=1$, as follows:

$$ \ddot{x}_1 + \omega_1^2 x_1 + \frac{1}{2} \omega_1^2 (3x_1^2 + x_2^2) + \omega_2^2 x_1 x_2 + \frac{1}{2} (\omega_1^2 + \omega_2^2) x_1 (x_1^2 + x_2^2) = 0, \quad (41) $$

$$ \ddot{x}_2 + \omega_2^2 x_2 + \frac{1}{2} \omega_2^2 (3x_2^2 + x_1^2) + \omega_1^2 x_1 x_2 + \frac{1}{2} (\omega_1^2 + \omega_2^2) x_2 (x_1^2 + x_2^2) = 0, \quad (42) $$
<figure><img src="image_11.png" /><figcaption>Fig. 8. Example 3—two dof nonlinear system with one mass.</figcaption></figure>

where $\omega_1^2 = k_1/m$, $\omega_2^2 = k_2/m$. The underlying linear system is uncoupled due to the orthogonal configuration of the two springs at rest.

The results presented here correspond to the following parameter values: $\omega_1 = 1$ and $\omega_2 = 2$, reflecting the presence of a two-order internal resonance, $\omega_2 = 2\omega_1$, of the underlying linear system at the static equilibrium point $X_0 = (0,0)$. The procedure outlined above was applied to generate the first two NNMs with an accurate ($m=4$)-Galerkin approximation for the selected amplitude ranges. Assumption H2 being not satisfied, the additive term $B(a, \Phi)$ is different from zero.

The behaviour of the first NNM is illustrated in Figs. 9–12. The second nonlinear mode is depicted in Figs. 13–16. We can observe that the influence of the even terms in the Fourier series increases with $a$. In both cases, the time histories of the periodic motions obtained using the proposed formulation (with the initial condition $\phi = 0$) and the direct numerical integration of Eqs. (41) and (42) are indistinguishable.

**5. Conclusion and further studies**

The NNMs have been formulated, as in the linear case, in terms of frequency and mode shape vector with the distinctive feature that these quantities are amplitude and total phase dependent. The formulation describes the NNM in terms of synchronous periodic oscillation as well as in terms of a two-dimensional invariant set of the dynamical system. So, the vibratory properties are distinguished from the geometrical properties.

The frequency and mode shape vector provide the solution to a nonlinear eigenvector–eigenvalue $2\pi$-periodic differential problem. As in the linear case, the eigenvalue (frequency) is always
<figure><img src="image_12.png" /><figcaption>Fig. 9. Example 3—behaviour of the first nonlinear normal mode. In plots (a)–(d): $m = 4$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Example 3—behaviour of the first nonlinear normal mode (continued). In plots (a) and (b): $m = 4$. In plot (b): $m = 0$ (continuous line), $m = 1, 3$ and $4$ (dashed lines).</figcaption></figure>
<figure><img src="image_14.png" /></figure>

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

<figure><img src="image_16.png" /><figcaption>Fig. 11. Example 3—behaviour of the first nonlinear normal mode (continued). In plot (a): invariant manifold in the phase subspace ($X_1, Y_1, X_2$) with $m = 4$. In plot (b): modal lines in the phase subspace ($X_1, Y_1$) ($m = 4$) for $a = 0.1$ and $0.41$. In plot (c): modal lines in the configuration space ($X_1, X_2$) ($m = 4$) for $a = 0.1$ and $0.41$.</figcaption></figure>

positive and can be expressed in terms of the eigenvector (mode shape). The eigenvector–eigenvalue $2\pi$-periodic problem is numerically solved using the classical Galerkin procedure in the space of $2\pi$-periodic vector functions. It should be noted that the zero-order solution reduced to the Szemplinska-Stupnicka [3] approach. Once the frequency is obtained, the period of the associated vibratory motion is easily deduced. It only depends on the amplitude. It seems that the presence of internal resonances in the underlying linear system does not require particular attention in these calculations.

The above formulation can be easily extended to damped autonomous nonlinear mechanical systems by introducing a scalar damping coefficient also depending upon the amplitude and the total phase. The computational Galerkin method must be modified accordingly [14].

It would be interesting for applications to express the general solution of system (1) under the form

$$X(t) = \sum_{l=1}^{n} \Psi_{l}(\mathbf{a}, \mathbf{\Phi}(t))a_{l} \cos \Phi_{l}(t),$$
<figure><img src="image_17.png" /><figcaption>Fig. 12. Example 3—time histories of the periodic motion for the first nonlinear normal mode: a = 0.44 and φ = 0. In plot (a): Φ versus time obtained by solving numerically Eq. (31). In plots (b) and (c): x₁ and x₂ versus time obtained from Eqs. (30) and (31) (dashed line) and by direct simulation of Eqs. (41) and (42) (continuous line).</figcaption></figure>

$$ \dot{\Phi}_l(t) = \Omega_l(a, \Phi_l(t)), \quad \Phi_l(0) = \phi_l, \quad \text{for } l = 1, \dots, n, $$

where the pair $(\Psi_l, \Omega_l)$ stands for the l-(coupled) mode, and $\mathbf{a} = (a_1, a_2, \dots, a_n)^T$, $\mathbf{\Phi} = (\Phi_1, \Phi_2, \dots, \Phi_n)$, denote the amplitude vector and the total phase vector, respectively. A promising approach has been proposed in Refs. [15,16] for solving nonlinear systems with wide-band random excitations.

## Appendix A

In this appendix, properties (17) and (18) of a well-defined periodic solution of Eqs. (14) and (16) are derived.

Multiplying Eq. (14) by $\Psi^T$ and using

$$ \Psi^T[M]\Psi = 1, \quad \Psi^T[M]\frac{\partial\Psi}{\partial\phi} = 0, \quad \Psi^T[M]\frac{\partial^2\Psi}{\partial\phi^2} = -\frac{\partial\Psi^T}{\partial\phi}[M]\frac{\partial\Psi}{\partial\phi}, $$
<figure><img src="image_18.png" /><figcaption>Fig. 13. Example 3—behaviour of the second nonlinear normal mode. In plots (a)–(d): $m = 4$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 14. Example 3—behaviour of the second nonlinear normal mode (continued). In plots (a) and (b): $m = 4$. In plot (b): $m = 0$ (continuous line), $m = 1$, 3 and 4 (dashed lines).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 15. Example 3—behaviour of the second nonlinear normal mode (continued). In plot (a): invariant manifold in the space subphase ($X_2, Y_2, X_1$) with $m = 4$. In plot (b): modal lines in the phase subspace ($X_2, Y_2$) ($m = 4$) with $a = 0.2$ and $0.8$. In plot (c): modal lines in the configuration space ($X_1, X_2$) ($m = 4$) with $a = 0.2$ and $0.8$.</figcaption></figure>

we obtain

$$ (\gamma^2(a, \phi) + 1)\Omega^2 \cos \phi + \frac{1}{2} \sin \phi \frac{\partial \Omega^2}{\partial \phi} = \frac{1}{a} \Psi^T F(\Psi a \cos \phi). \quad (43) $$

Let $\Omega^2(a, \phi^*)$ be some extremum value of $\Omega^2(a, \cdot)$ on $[0, 2\pi]$. Multiplying Eq. (43) by $\cos(\phi^*)$ gives

$$ (\gamma^2(a, \phi^*) + 1)\cos^2 \phi^* \Omega^2 = \frac{1}{a^2} \Psi^T a \cos \phi^* F(\Psi a \cos \phi^*), \quad (44) $$

where we have used that $(\partial\Omega^2/\partial\phi)(a, \phi^*) = 0$. By H1, the right-hand side of Eq. (44) is non-negative. Consequently, any extremum value will be also non-negative and the function $\Omega^2(a, \cdot)$ will be non-negative everywhere on $[0, 2\pi]$.

Setting $Y = \Omega^2 \sin^2 \phi$, Eq. (43) becomes

$$ \beta(a, \phi) Y + \frac{\partial Y}{\partial \phi} = \frac{2}{a} \sin \phi \Psi^T F(\Psi a \cos \phi), \quad (45) $$
<figure><img src="image_21.png" /><figcaption>Fig. 16. Example 3—time histories of the periodic motion for the second NNM: a = 0.8 and φ = 0. In plot (a): Φ versus time obtained by solving numerically Eq. (31). In plots (b) and (c): x₁ and x₂ versus time obtained from Eqs. (30) and (31) (dashed line) and by direct simulation of Eqs. (41) and (42) (continuous line).</figcaption></figure>

and the constant variation formula gives

$$Y(\phi) = \frac{2}{a} \left( \exp\left(-\int_{0}^{\phi} \beta(a, v) \, dv\right) \right) \left( \int_{0}^{\phi} \exp\left(\int_{0}^{\sigma} \beta(a, v) \, dv\right) I(a, \sigma) \, d\sigma \right). \quad (46)$$

Owing to the symmetry properties of a well-defined periodic solution for $\Psi(\Psi(a, \phi) = \Psi(a, -\phi) = \Psi(a, \phi + \pi))$, Eq. (46) makes sense as a $\pi$-periodic solution of Eq. (45) such that $Y(\pi) = Y(0) = 0$. $I(a, \phi)$ (as well as $\beta(a, \phi)$) is in fact a regular $\pi$-periodic function with respect to $\phi$, having a Fourier series of the form $I(a, \phi) = \sum_{k=1}^{\infty} I_{2k}(a) \sin 2k\phi$, so that $\int_{0}^{\pi} \exp(\int_{0}^{\sigma} \beta(a, v) \, dv) I(a, \sigma) \, d\sigma = 0$. Finally, expression (18) makes sense for $\phi = 0$ and $\pi$, and can be expanded into $\Omega^2(a, \phi) = \sum_{k=0}^{\infty} \Omega_{2k}^2 \cos 2k\phi$.

References

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

[2] R.M. Rosenberg, On nonlinear vibrations of systems with many degree of freedom, *Advances in Applied Mechanics* 242 (9) (1966) 155–242.
[3] W. Szemplinska-Stupnicka, in: *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-448.

[5] C. Touze, A normal form approach for nonlinear normal modes, in: S. Bellizzi, B. Cochelin (Eds.), *Quelques aspects des modes non-linéaires en vibrations*, Publications du LMA 156, 2003, pp. 15-38, ISSN 0750-7356.

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

[7] 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.

[8] B. Cochelin, Modes non-linéaires de structures élastiques: quelques définitions et méthodes de calcul, in: S. Bellizzi, B. Cochelin (Eds.), *Quelques aspects des modes non-linéaires en vibrations*, Publications du LMA, 156, 2003, pp. 1-14, ISSN 0750-7356.

[9] M. Urabe, Galerkin's procedure for nonlinear periodic systems, *Archs. Ration. Mech. Analysis* 20 (1965) 120-152.

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

[11] L. Azrar, R. Benamar, R.G. White, A semi-analytical approach to the nonlinear dynamic response problem of S-S and C-C beams at large vibration amplitudes part 1: general theory and application to the single mode approach to free and forced vibration analysis, *Journal of Sound and Vibration* 224 (2) (1999) 183-207.

[12] P. Ribiero, M. Petyt, Nonlinear vibration of plates by the hierarchical finite element and continuation methods, *International Journal of Mechanical Sciences* 41 (1999) 437-459.

[13] C. Touze, O. Thomas, Reduced order modeling for a cantilever beam subjected to harmonic forcing, 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. 165-168. ISBN: 2-86 834-119-5.

[14] S. Bellizzi, R. Bouc, A new formulation for the existence and calculation of nonlinear modes, 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. 25-28, ISBN: 2-86 834-119-5.

[15] S. Bellizzi, R. Bouc, Analysis of multi-degree of freedom strongly nonlinear mechanical systems with random input. Part I: nonlinear modes and stochastic averaging, *Probabilistic Engineering Mechanics* 14 (1999) 229-244.

[16] S. Bellizzi, R. Bouc, Analysis of multi-degree of freedom strongly nonlinear mechanical systems with random input. Part II: equivalent linear system with random matrices and power spectral density matrix, *Probabilistic Engineering Mechanics* 14 (1999) 245-256.