# Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes

Rémi Arquier, Sergio Bellizzi, Robert Bouc, Bruno Cochelin

► To cite this version:

Rémi Arquier, Sergio Bellizzi, Robert Bouc, Bruno Cochelin. Two methods for the computation of nonlinear modes of vibrating systems at large amplitudes. Computers & Structures, 2006, 84 (24-25), pp.1565-1576. 10.1016/j.compstruc.2006.01.011. hal-01580935

HAL Id: hal-01580935

https://hal.science/hal-01580935

Submitted on 3 Sep 2017

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 Attribution 4.0 International License
Two methods for the computation of nonlinear modes
of vibrating systems at large amplitudes

Rémi Arquier, Sergio Bellizzi, Robert Bouc, Bruno Cochelin

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

The aim of this paper is to present two methods for the calculation of the nonlinear normal modes of vibration for undamped non-
linear mechanical systems: the time integration periodic orbit method and the modal representation method. In the periodic orbit
method, the nonlinear normal mode is obtained by making the continuation of branches of periodic orbits of the equation of motion.
The terms “periodic orbits” means a closed trajectory in the phase space, which is obtained by time integration. In the modal represen-
tation method, the nonlinear normal mode is constructed in terms of amplitude, phase, mode shape, and frequency, with the distinctive
feature that the last two quantities are amplitude and total phase dependent. The methods are compared on two DOF strongly nonlinear
systems.

Keywords: Vibration; Normal mode; Nonlinear structure; Periodic motion

**1. Introduction**

Extending the concept of normal modes of vibration to
the case where the restoring forces contain nonlinear terms,
has been a challenge to many authors. This has led to the
so-called nonlinear normal modes (NNMs) which have
great potential for applications in nonlinear vibrating sys-
tems. For instance, a damped system exited by harmonic
forcing will have its resonances close to the NNMs, and
this is a first obvious motivation for their computation. It
is also now established that the knowledge of the NNMs,
together with their bifurcations, can be very helpful to
understand the dynamics of a nonlinear system [1]. Some
important phenomenon such as the localisation of the
motion [2], the interaction between modes [3], the pumping
of energy of a linear system by a pure nonlinear one [4],
can be nicely explained using the NNM concept. The
NNMs are also important for the modal controllability

of nonlinear systems [5]. Finally, eventhough the principle
of superposition does not hold for nonlinear equations, the
NNM can be useful to generate effective reduced-order
models for multi-degree-of-freedom nonlinear systems
[6,7].

Following the pioneer work by Rosenberg [8] on conser-
vative systems, several attempts have been made to develop
methods for the calculation of nonlinear normal modes.
Without entering into details, we mention here several clas-
ses of techniques which have aimed this goal. They are the
harmonic balance approach [9–12], the normal form theory
[13–15], many perturbation techniques [16,17] such as the
famous multiple scale analysis, and the invariant manifold
method [18,19] which led to a new definition of NNMs,
extending the concept to nonconservative systems. Many
of the above mentioned methods are based on some kind
of perturbation expansions with a truncature of the series
after the first few leading order terms. They have the
advantage to provide analytical expressions of the NNMs,
but on the other hand, the drawback to be limited to weak
nonlinearities or small amplitudes. Is it now evident, see
[7,20] for instance, that numerical methods should be developed to complete the analytical ones. These numerical methods should be able to explore the NMN at large amplitude of vibration, and possibly, to detect all kind of bifurcations along a NNM. The aim of this paper is to present and compare two approaches that fall in that category. The first one is the time integration of a periodic orbit and the second, the modal representation method. In both of them, we compute a one dimensional family of periodic solutions that are parametrized either by the level of mechanical energy or by the amplitude of the orbits. These families of periodic orbits provide a two dimensional invariant subspace of the phase space which has the property to pass through the equilibrium point (zero amplitude orbits) and to be tangent to the linear mode of the underlying linear system (small amplitude orbits). Accordingly, these families of orbits correspond to the NNM as defined by Shaw and Pierre [18,19].

The paper is organized as follows: the mechanical framework is given in Section 2. Section 3 is devoted to the first method which is the time integration of periodic orbit. The NNMs are computed here by making the continuation of branches of periodic orbits that are parametrized by the energy level. The closed orbits are computed using an “exact energy conserving” time integration algorithm [21]. This transforms the boundary value problem into an algebraic one which depends on a free parameter, and which is solved with a semi-analytical continuation technique: the so-called asymptotic-numerical method (ANM) [22]. This approach provides the period of the oscillations on the NNMs, and a time discrete representation of the orbits. The determination of the stability and of the bifurcation of the NNM is well established with this numerically oriented approach [23]. It will presented in a forthcoming paper. Section 4 is devoted to the second method called “the modal representation method”. As in the linear case, an expression is developed for the NNM 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 dynamics 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π-periodic nonlinear eigenvalue problem from which a numerical Galerkin procedure is developed for approximating the NNMs. This formulation allows us to characterize the similar NNMs. It leads also to an analytical (parametric) expression of the invariant manifold and it permits to compute the NNM even in the case of resonance relations between the eigenvalues of the linearized system. The extension of the formulation in the case of damped autonomous mechanical systems is considered in [24]. Finally, in Section 5, these two methods are compared on a benchmark problem with two Green-Lagrange springs, which is representative of geometrically nonlinear thin structures such as plates and shells.

## 2. Mechanical framework

In this study, we consider the undamped autonomous nonlinear n degrees of freedom (n-DOF) mechanical system

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

where M is the mass matrix, and F(U) is the vector of restoring forces including linear and nonlinear terms. The overdots stand for temporal derivatives. The following assumptions will be made throughout this study:

* H1: M is a symmetrical positive definite matrix;

* H2: $F(U) = \frac{\partial W}{\partial U}(U)$ where $W(U)$ is a scalar potential energy function of the vector $U \in R^n$, $W$ has a continuous second derivatives, is positive, and admits a local minimum at $U = 0$. (As usual $\frac{\partial W}{\partial U}(U)$ will be denoted $W_{,U}(U)$).

The conservative system (1) has a first integral corresponding to the conservation of the total energy $E$, which is the sum of the kinetic energy and the potential energy $W$ i.e.

$$ E(U(t)) = \frac{1}{2} \dot{U}(t)^{\mathrm{T}} M \dot{U}(t) + W(U(t)). \quad (2) $$

The linear equation

$$ M\ddot{U}(t) + [W_{,UU}(0)]U(t) = 0, \quad (3) $$

where $[W_{,UU}(0)]$ denotes the Hessian matrix for the function $W$ at $U = 0$, will be called the underlying linear system (or linearized system) associated with the nonlinear equation (1).

It should be noted that the framework (1) includes the equations of motion of elastic thin structures with geometrical nonlinearities such as shells, plates, beams and cables. The continuous model should be discretised using a classical Ritz or finite element method.

In the following, we focus on the periodic solution of (1). If there is no internal resonance (the eigenfrequencies of the linearized system (3) are no commensurable), the conservative system (1) possesses at least n two dimensional families of periodic solution around the stable origin $U = 0$. These two dimensional families of periodic orbits allow description of two dimensional invariant manifolds of the phase space, corresponding to the NNM, as defined by Shaw and Pierre [18].

## 3. Periodic orbits method

The numerical computation of periodic orbits has already been addressed in textbooks [23,25] either for the calculation of isolated orbits or for the continuation of a family of orbits. The most popular method is the so-called shooting method which consists in finding a suitable initial condition, that induces a closed trajectory in the phase space. This leads to a boundary value problem [23] where
the boundary condition is in fact a periodicity condition. Since the physical system is nonlinear, iterative process like Newton-Raphson method is intensively used to find the initial condition. It should be noted that each Newton-Raphson iteration requires the entire computation of the orbit by direct time integration. The whole process is thus relatively heavy and can fail when the periodic solutions are unstable.

Shooting techniques have already been used to determine the forced response of nonlinear damped systems [26]. Once one orbit has been computed for a given value of the forcing term frequency, a pseudo-arc length continuation procedure is performed to obtain a diagram response of the amplitude response with respect to the forcing frequency.

Our method keeps the essential aspects of continuation of periodic orbits, but with some adaptation for application to the NNM of conservative systems at large amplitudes.

* We do not consider any forcing terms nor damping in the governing equation. Hence, the system has a first integral corresponding to the conservation of the total energy. This brings some additional troubles in the numerical determination of the periodic orbits, that are addressed in Section 3.1.

* As an alternative to shooting where only the initial point $U(0), \dot{U}(0)$ is unknown, we prefer a global approach where the unknown is the values of the vector $U(t), \dot{U}(t)$ at the temporal grid point. This global approach is known to have a better convergence and is more robust than the shooting one [27]. It is presented in Section 3.2.

* For the continuation of the closed orbits, we use the ANM which is more convenient than the classical Newton-Raphson for following branches and for detecting bifurcations [28].

For convenience, Eq. (1) is re-written as a first-order differential equation

$$ \mathbf{R}(\mathbf{Z}) = \mathbf{A}\dot{\mathbf{Z}} + \mathbf{f}(\mathbf{Z}) = \mathbf{0}, \quad (4) $$

where $\mathbf{Z} = [\mathbf{U}^T, \mathbf{V}^T]^T$ is the 2n state vector, $\mathbf{V} = \dot{\mathbf{U}}$ denotes the velocity,

$$ \mathbf{A} = \begin{bmatrix} \mathbf{I} \\ \mathbf{M} \end{bmatrix} \quad \text{and} \quad \mathbf{f}(\mathbf{Z}) = \begin{bmatrix} -\mathbf{V} \\ \mathbf{F}(\mathbf{U}) \end{bmatrix}. $$

### 3.1. Periodic orbits for a conservative system

The terms “periodic solution” or “periodic orbit” mean a closed trajectory in the phase space. It is defined as follows:

$\mathbf{Z}(t)$ is a periodic orbit of Eq. (4), with the period $T$, if it solves the system

$$ \mathbf{S}_{\text{under}}(\mathbf{Z}(t), T) = \begin{cases} \mathbf{A}\mathbf{Z}(t) + \mathbf{f}(\mathbf{Z}(t)) = \mathbf{0} & \text{for } 0 < t < T, \\ \mathbf{Z}(0) - \mathbf{Z}(T) = \mathbf{0}. & \end{cases} \quad (5) $$

Eq. (5) is the governing equation that must be satisfied over the time interval $0 < t < T$, where $T$ is the unknown period of the periodic orbit. The second equation is the periodicity condition which forces the first state $\mathbf{Z}(0)$ to be equal to the last $\mathbf{Z}(T)$ to close the trajectory.

For autonomous conservatives systems with a first integral, Eq. (5) defines a two dimensional family of periodic orbits. For the numerical solution of Eq. (5) it is mandatory to add other conditions in order to define unique closed orbits belonging to these two dimensional spaces. Otherwise, numerical problem such as noninvertible matrix occurs.

* The first condition is the classical “phase condition” [23,25] which forces the initial point $\mathbf{Z}(0)$ to be unique for a given periodic orbit. Here, the component $i$ of the initial state vector $\mathbf{Z}_0$ is set to zero. Usually, this component corresponds to a velocity degree of freedom, so we have $i \in [n+1, 2n]$,

$$ \{\mathbf{Z}(0)\}_i = 0. \quad (6) $$

* The second condition, named “Energy condition”, prescribes the value of the total mechanical energy on the periodic orbit. By this way the value $E_0$ “select” one orbit on the NNM i.e.

$$ E(\mathbf{Z}(0)) = E_0. \quad (7) $$

The new system, denoted $\mathbf{S}_{\text{over}}$, obtained by grouping Eqs. (5)-(7):

$$ \mathbf{S}_{\text{over}}(\mathbf{Z}(t), T, E_0) = \begin{cases} \mathbf{A}\mathbf{Z}(t) + \mathbf{f}(\mathbf{Z}(t)) = \mathbf{0} & \text{for } 0 < t < T, \\ \mathbf{Z}(0) - \mathbf{Z}(T) = \mathbf{0}, \\ \{\mathbf{Z}(0)\}_i = 0, \\ E(\mathbf{Z}(0)) = E_0 \end{cases} \quad (8) $$

defines one isolated orbit<fnref n="1" /> of Eq. (4). It is possible to solve numerically $\mathbf{S}_{\text{over}}$ (after application of a suitable discretisation scheme), and obtain periodic orbits for different values of $E_0$. However, even though the solution of $\mathbf{S}_{\text{over}}$ is unique, the system $\mathbf{S}_{\text{over}}$ is overdetermined, there is one more equation than unknowns. This type of system could be numerically solved by Least Square Methods but this situation is not suitable for the continuation method we want to use, and not convenient when the number of unknown is large.

To remove this drawback we follow a solution proposed in [29] and extended in [30] for systems with multiple first integrals. It consists in replacing the conservative system (1) by a nonconservative system which contains the

<footnote n="1">¹ One isolated orbit if we are on a regular point on the branch.</footnote>
dissipative term $A \frac{\partial E}{\partial Z}(Z(t))$ into the governing equation. This leads to the following new system, denoted $S_{\text{square}}$,

$$ S_{\text{square}}(\mathbf{Z}(t), T, \beta, E_0) = \left\{ \begin{array}{l} \mathbf{A}\dot{\mathbf{Z}}(t) + \mathbf{f}(\mathbf{Z}(t)) + \beta\mathbf{A}\frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t)) = 0 \quad \text{for } 0 < t < T, \\ \mathbf{Z}(0) - \mathbf{Z}(T) = 0, \\ \{\mathbf{Z}(0)\}_i = 0, \\ E(\mathbf{Z}(0)) = E_0. \end{array} \right. \quad (9) $$

with a new scalar unknown $\beta$ which controls the amount of dissipation. The dissipative term is the gradient of the total energy $E$. It is acting like a force which is always perpendicular to the orbit, and which tends to bring $Z$ toward the origin. The new system equation (9) can have a periodic orbit only if $\beta = 0$. Indeed, if $\mathbf{Z}(0) - \mathbf{Z}(T) = 0$, we have

$$ \begin{aligned} E(\mathbf{Z}(T)) - E(\mathbf{Z}(0)) &= \int_0^T \frac{dE}{dt}(\mathbf{Z}(t)) dt = \int_0^T \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t))^T \dot{\mathbf{Z}}(t) dt \\ &= -\int_0^T \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t))^T (\mathbf{A}^{-1}\mathbf{f}(\mathbf{Z}(t)) + \beta \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t))) dt \\ &= -\int_0^T \left( \begin{bmatrix} -\mathbf{F}(\mathbf{U}) \\ \mathbf{M}\mathbf{V} \end{bmatrix}^T \begin{bmatrix} I \\ \mathbf{M}^{-1} \end{bmatrix} \begin{bmatrix} -\mathbf{V} \\ \mathbf{F}(\{\mathbf{U}\}) \end{bmatrix} \right) + \beta \left\| \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t)) \right\|^2 dt \\ &= -\int_0^T \mathbf{F}(\mathbf{U})^T \mathbf{V} - \mathbf{V}\mathbf{F}(\mathbf{U}) + \beta \left\| \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t)) \right\|^2 dt \\ &= -\beta \int_0^T \left\| \frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t)) \right\|^2 dt = 0, \end{aligned} \quad (10) $$

which proves that $\beta = 0$ whatever the value of $\frac{\partial E}{\partial \mathbf{Z}}(\mathbf{Z}(t))$. Hence, the solutions given by Eq. (9) correspond to the solutions of the original conservative system (8). The advantage of Eq. (9) over (8) is that the Jacobian matrix is square and non-singular.<fnref n="2" />

## 3.2. Time discretisation and global solving

The implicit Newmark method is the most classical numerical schemes for the time integration in structural dynamics. In the nonlinear regime, it has been shown that this scheme can be unstable for conservative systems [31]. Simo and Tarnow have proposed an energy conserving discretisation scheme to remove this drawback [21]. They introduce a subtle modification of the internal forces in the discrete equations, which induce a conservation of the total energy between two time steps. Such a conservative time stepping scheme is well adapted to find closed orbits in a system which has no forcing, nor dissipation. Indeed, in presence of numerical damping for instance, the trajectory could not come back at the same energy level after a period of integration, unless very small steps are used in order to limit the numerical dissipation. In regards of these considerations, we have chosen the Simo scheme

to re-write the system (9) in a “time discretised manner”. We shall not enter in the details of the Simo integration scheme, and we will now describe how we obtain one periodic orbit.

According to our numerical experience, see also [27], we preferred to use a global solution technique instead of the more classical shooting method. By global solution technique, we mean that the unknown vector $\mathbf{Z}_{op}$ contains all the state vectors $\mathbf{Z}_{k,k=0,...,m}$ of the system for each time steps $t_{k,k=0,...,m}: \mathbf{Z}_{op} = [\mathbf{Z}_0, \mathbf{Z}_1, \mathbf{Z}_2, ..., \mathbf{Z}_m]^T$. With $\mathbf{Z}_k = \mathbf{Z}(t = t_k)$. Hence, if we have $m + 1$ time steps for one orbit, we write the governing equation $m$ time in the same system. This leads to a large and very sparse system, which can be written as follows:

$$ S_{\text{square}}(\mathbf{Z}_{\text{op}}, T, \beta, E_0) = \left\{ \begin{array}{ll} \mathbf{R}_d(\mathbf{Z}_0, \mathbf{Z}_1, T) + \beta \mathbf{D}_d(\mathbf{Z}_0, \mathbf{Z}_1, T) & = 0, \\ \mathbf{R}_d(\mathbf{Z}_1, \mathbf{Z}_2, T) + \beta \mathbf{D}_d(\mathbf{Z}_1, \mathbf{Z}_2, T) & = 0, \\ \vdots & \vdots \\ \mathbf{R}_d(\mathbf{Z}_{m-1}, \mathbf{Z}_m, T) + \beta \mathbf{D}_d(\mathbf{Z}_{m-1}, \mathbf{Z}_m, T) & = 0, \\ \mathbf{Z}_0 - \mathbf{Z}_m & = 0, \\ \{\mathbf{Z}_0\}_i & = 0, \\ E_d(\mathbf{Z}_0, \mathbf{Z}_1, T) & = E_0, \end{array} \right. \quad (11) $$

where $\mathbf{R}_d$ is the discretised form of the behaviour equation and $\mathbf{D}_d$ is discretised form of the dissipative term. Note also that we use the classical time normalisation procedure in order to have a constant number of time steps. This normalisation consists in replacing the real time step $dt$ by the following expression:

$$ dt = T/m, $$

where $T$ is the period and $m$ the number of time steps which remains constant. This explains why $\mathbf{R}_d$, $\mathbf{D}_d$ and $E_d$ depend on the period $T$.

Solving this system for a given value of $E_0$ permits to obtain one periodic orbit on the NNM. It is possible to do this with an initial good estimate and some Newton corrections.

## 3.3. Numerical construction of the NNM

For the construction of one NNM, we start by finding a periodic orbit near the linear domain where nonlinear terms are negligible: we use the linear mode as a first estimate, then we apply a classical correction method to return on the solution branch. By this way we obtain the first point (one orbit) for our continuation. Next, we follow the branch by the ANM.

The ANM is a continuation method to solve nonlinear algebraic problems $S(\mathbf{U}) = 0$, $\mathbf{S} \in \mathbb{R}^n$ and $\mathbf{U} \in \mathbb{R}^{n+1}$ [22, 32]. This method is based on power series expansions of the unknowns $\mathbf{U}(s)$ with respect to the path parameter $s$ and it consists in generating a succession of branches, instead of a sequence of points.

<footnote n="2">² Except at bifurcation points of course.</footnote>
Assuming $\mathbf{U} = [\mathbf{Z}_{op}, T, \beta, E_0]^T$ the system (11) can be re-written in a quadratic form<fnref n="3" /> like

$$S(\mathbf{U}) = \mathbf{L}(\mathbf{U}) + \mathbf{Q}(\mathbf{U}, \mathbf{U}), \qquad (12)$$

where $\mathbf{L}$ is linear, and $\mathbf{Q}$ quadratic. Denoting $\mathbf{U}_0$ a known solution point, introducing the following expansion

$$\mathbf{U}(s) = \mathbf{U}_0 + \sum_{i=1}^{n} s^i \mathbf{U}_i, \qquad (13)$$

in Eq. (12) and equating the terms at the same power of $s$, results in the decomposition of the nonlinear system in a sequence of linear systems that will be solved successively. To complete the system, we need to add a condition which is the definition of $a$, the path parameter. For instance we can use the classical pseudo-arc length parameter $s$ defined by

$$\mathbf{U}^T \mathbf{U}_1 = s. \qquad (14)$$

Hence, the linear problems can be written as follows:

$$\text{order 1:} \quad \begin{cases} \mathbf{L}_t(\mathbf{U}_1) = 0, \\ \mathbf{U}_1^T \mathbf{U}_1 = 1, \end{cases} \qquad (15)$$

$$\text{order } p: \quad \begin{cases} \mathbf{L}_t(\mathbf{U}_p) = \mathbf{F}_p^{nl}, \\ \mathbf{U}_p^T \mathbf{U}_1 = 0, \end{cases} \qquad (16)$$

where $\mathbf{L}_t(\cdot)$ is a tangent operator at $\mathbf{U} = \mathbf{U}_0: \mathbf{L}_t(\cdot) = \mathbf{L}(\cdot) + \mathbf{Q}(\cdot, \mathbf{U}_0) + \mathbf{Q}(\mathbf{U}_0, \cdot)$.

Each resolution of Eq. (16) at order $p$ give the terms $\mathbf{U}_p$. The domain of validity of the solution is defined a posteri-
ori by analysing the convergence radius of the series, and
so, a part of a branch is calculated for each step. Advanta-
ges of this method mainly lie in three points: analytical
solutions are obtained, only one matrix inversion is needed
for each part of branch; eventually, robust continuation is
performed, using optimal steps, and this method is much
more reliable than classical incremental-iterative methods
such as Newton-Raphson. Bifurcation points [28] are also
easier to detect.

In summary, the ANM provides a succession of contin-
uous representation of the family of periodic orbits. By this
way, we obtain the two dimensional invariant surface of
the state space, introduced by Shaw and Pierre [18,19].

## 4. Modal representation method

Here also the NNM will be defined in the framework of periodic solutions of Eq. (1) and the following properties of periodic functions will be used. If $\mathbf{U}^*(t)$ is a periodic solution with period $T > 0$ of Eq. (1) (i.e. $\mathbf{U}^*(t) = \mathbf{U}(t+T)$, $\forall t$), then the function $\mathbf{U}^{**}$ defined by $\mathbf{U}^{**}(t) = \mathbf{U}^*(-t)$ is also a $T$-periodic solution for Eq. (1). Hence this section will be devoted to the $T$-periodic solution of Eq. (1) with

symmetry $\mathbf{U}^*(t) = \mathbf{U}^*(-t)$. This property implies that the associated Fourier series will contain only cosine terms.

### 4.1. Linear case

If $\mathbf{F}(\mathbf{U}) = \mathbf{KU}$ where $\mathbf{K}$ is a symmetrical positive definite matrix, the general solution of Eq. (1) can be expressed as

$$\mathbf{U}(t) = \sum_{p=1}^{n} a_p \psi_p \cos(\omega_p t + \phi_p),$$

where $(\omega_p, \psi_p)$ solve the eigenproblem

$$\mathbf{K}\psi_p = \mathbf{M}\psi_p\omega_p^2 \text{ with } \psi_p^T\mathbf{M}\psi_p = 1 \text{ and } \psi_p^T\mathbf{M}\psi_l = 0, \quad p \neq l. \qquad (17)$$

The contribution of a mode is given by (where the subscript $p$ has been omitted)

$$\mathbf{U}(t) = a\psi \cos(\omega t + \varphi)$$

and can be re-written as

$$\mathbf{U}(t) = a\psi \cos(\Phi(t)) \qquad (18)$$

with

$$\begin{cases} \dot{\Phi}(t) = \omega, \\ \Phi(0) = \varphi. \end{cases} \qquad (19)$$

The scalar function $\Phi(\cdot)$ (or total phase) describes the time evolution on the modal line (a straight line) in the configuration space.

### 4.2. Nonlinear case with odd restoring force

In addition to H1 and H2, we assume throughout this
subsection that

$$\bullet \text{ H3: } \mathbf{F}(\mathbf{U}) = -\mathbf{F}(-\mathbf{U}).$$

As suggested by the linear case, a normal mode is defined as a periodic solution to Eq. (1) and is sought in the form

$$\mathbf{U}(t) = a\mathbf{\Psi}(a, \Phi(t)) \cos(\Phi(t)) \qquad (20)$$

with

$$\begin{cases} \dot{\Phi}(t) = \Omega(a, \Phi(t)), \\ \Phi(0) = \varphi. \end{cases} \qquad (21)$$

The vector $\mathbf{\Psi}$ and the scalar $\Omega$ are now functions which depend on two variables and they are $2\pi$-periodic with respect to the second variable. By analogy with the linear case, the functions, $\mathbf{\Psi}$ and $\Omega$, will be referred to here as the modal vector (function) and the resonance frequency (function) of the nonlinear normal mode, respectively. The one dimensional differential equation (21) governing the total phase motion $\Phi$ will define the dynamics of the periodic response and the scalar quantities $a$ ($a > 0$) and $\varphi$ ($\varphi \in [0, 2\pi]$) will set the initial conditions of the vibration in mode motion on the phase space. To ensure that the parameter

<footnote n="3">³ The details of this re-written form will be described in a forthcoming paper.</footnote>
a appropriately characterizes the amplitude of the vibration in mode motion, a normalization condition on the modal vector **Ψ** is required. In this study (and without loss of generality), we will adopt the following condition:

$$ \mathbf{\Psi}^T(a, \phi)\mathbf{M}\mathbf{\Psi}(a, \phi) = 1 \quad \forall (a, \phi) \in \mathbb{R}^+ \times [0, 2\pi]. \quad (22) $$

Differentiating Eq. (20), substituting into Eq. (1) and eliminating the time variable in the resulting equation, Eq. (1) reduces to the following differential equation:

$$ \begin{aligned} \mathbf{L}(\Omega^2, \mathbf{\Psi}; \phi) &+ \frac{1}{a}\mathbf{F}(\mathbf{\Psi}a \cos \phi) \\ &= \mathbf{M}\mathbf{\Psi}(\Omega^2 \cos \phi + \frac{1}{2}(\Omega^2)_{,\phi} \sin \phi), \end{aligned} \quad (23) $$

where the differential operator **L** is given by

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

and $\phi, \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. (23), which can be viewed as an extension of the eigenvalue–eigenvector problem (17), will be named nonlinear eigenvalue–eigenvector problem and will be used together with the normalization condition (22) to characterize the nonlinear normal mode of motion, i.e. to characterize the modal vector **Ψ** and the square frequency $\Omega^2$. Thus, the modal vector **Ψ** and the frequency $\Omega$ will be viewed here as amplitude and total phase dependent and will be searched for as $2\pi$-periodic functions with respect to $\phi$ for fixed a. More specifically, according to the symmetry assumption H3, these functions will be searched for as even periodic functions with respect to $\phi$ with period $\pi$. It follows that the Fourier series of $\mathbf{U}(t)$ in terms of $\mathbf{\Phi}(t)$ will contain only odd cosines terms.

For fixed $a \in \mathbb{R}^+$ and under H1–H3, it can be shown (see [33]) that there exist $n$ $2\pi$-periodic functions, $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$ with $\Omega_p^2(a, \cdot) > 0$, which solve Eqs. (23) and (22) with

(1) each solution $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$ is unique in some neighbourhood of $a=0$ and $(\omega_p^2, \psi_p)$ where the pair $(\omega_p^2, \psi_p)$ characterizes a normal mode of the underlying linear system;

(2) the period of the vibration in mode motion $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$ depends only on the amplitude $a$ and it is given by

$$ T_p(a) = 2 \int_0^\pi \frac{1}{\Omega_p(a, \phi)} d\phi; \quad (25) $$

(3) as $a \to 0$, the nonlinear normal mode $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$ tends toward the linear one $(\omega_p, \psi_p)$.

For each pair $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$, Eq. (20) together with $\dot{\mathbf{U}}(t) = a\Omega(a, \mathbf{\Phi}(t))(\mathbf{\Psi}_{,\phi}(a, \mathbf{\Phi}(t))\cos \mathbf{\Phi}(t)) - \mathbf{\Psi}(a, \mathbf{\Phi}(t))\sin \mathbf{\Phi}(t)$, \quad (26)

define a “synchronous” periodic oscillation [8] (nonlinear normal mode). 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 NNM in the framework of invariant manifold [18–20] in the phase space. For each pair $(\Omega_p^2(a, \cdot), \mathbf{\Psi}_p(a, \cdot))$, the invariant manifold coincides with the set, in the phase space, defined by Eqs. (20) and (26) taking the initial conditions $a$ and $\phi$ to be independent variables.

An illustrative example: Let us consider the following 2 d.o.f. nonlinear system governed by the equations of motion

$$ \ddot{u}_1 + (k+1)u_1 - u_2 + r(u_1 - u_2)^2(u_1 - u_2) = 0, \quad (27) $$

$$ \ddot{u}_2 - u_1 + (k+1)u_2 + r(u_1 - u_2)^2(u_2 - u_1) = 0, \quad (28) $$

where $k$ and $r$ are positive constant characterizing the linear and nonlinear stiffnesses, respectively.

It is easy to show that the following two pairs:

$$ \begin{cases} \Omega_1^2(a, \phi) = k, \\ \mathbf{\Psi}_1(a, \phi) = (1/\sqrt{2}, 1/\sqrt{2})^T, \end{cases} \quad \begin{cases} \Omega_2^2(a, \phi) = k + 2 + ra^2(3 + \cos 2\phi), \\ \mathbf{\Psi}_2(a, \phi) = (1/\sqrt{2}, -1/\sqrt{2})^T \end{cases} $$

solve the associated eigenvalue–eigenvector problem (23) and (22). This two pairs define two nonlinear normal modes of the system (27) and (28). Each nonlinear normal mode can be viewed as the continuation of the linear normal modes of the underlying linear system. The first mode is not affected by the nonlinearity whereas the nonlinearity only affects the dynamics $(\mathbf{\Psi}_2(a, \phi)$ is a constant vector) of the second mode. Both modes are similar.

### 4.3. General case

When $\mathbf{F}(\mathbf{U}) \neq -\mathbf{F}(-\mathbf{U})$, the periodic solution (to characterize the nonlinear normal modes) are sought in the form

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

with

$$ \begin{cases} \dot{\mathbf{\Phi}}(t) = \mathbf{\Omega}(a, \mathbf{\Phi}(t)), \\ \mathbf{\Phi}(0) = \phi, \end{cases} \quad (30) $$

where the vectors **B** and **Ψ**, and the scalar $\Omega$ are functions which depend on two variables and are $2\pi$-periodic with respect to the second variable. The term depending on **B** has been added in Eq. (29) to balance the even cosine terms which can appear in the Fourier series of $\mathbf{U}(t)$ in terms of $\mathbf{\Phi}(t)$ when assumption H3 is not satisfied. The function **B** will be referred to here as the *bias term* of the nonlinear normal mode. Hence **B** will be viewed (like **Ψ** and $\Omega$) as amplitude and total phase dependent, and will be searched for an even periodic function with respect to $\phi$ for fixed a.

Differentiating Eq. (29), substituting into Eq. (1) and eliminating the time variable in the resulting equation, Eq. (1) reduces now to the following differential equation:
$$ 
\mathbf{L}(\Omega^2, \mathbf{\Psi}; \phi) + \frac{1}{a}\mathbf{F}(\mathbf{\Psi}a \cos \phi + a\mathbf{B}) + \Omega^2 \mathbf{M}\mathbf{B}_{,\phi\phi} \\ 
+ \frac{1}{2}(\Omega^2)_{,\phi} \mathbf{M}\mathbf{B}_{,\phi} = \mathbf{M}\mathbf{\Psi}(\Omega^2 \cos \phi + \frac{1}{2}(\Omega^2)_{,\phi} \sin \phi), \quad (31) $$

where the differential operator $L$ is still given by Eq. (24). Here also 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. (31), which is an extension of the eigenvalue–eigenvector problem (23) will be used together with the normalization condition (22) to characterize the nonlinear normal mode of motion i.e. to characterize the modal vector $\mathbf{\Psi}$, the square frequency $\Omega^2$ and the bias term $\mathbf{B}$.

Apparently, we have to solve, for fixed $a$, two Eqs. (31) and (22) for three unknowns $(\Omega^2(a, \cdot), \mathbf{\Psi}(a, \cdot), \mathbf{B}(a, \cdot))$. Recalling the periodicity properties of the functions $\Omega^2$ and $\mathbf{\Psi}$ (as established in the previous section), the bias term $\mathbf{B}$ has to be a $2\pi$-periodic function with only even cosine terms and we can only look for even $\pi$-periodic solutions $(\Omega^2(a, \cdot), \mathbf{\Psi}(a, \cdot), \mathbf{B}(a, \cdot))$. Consequently using the decomposition

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

where $\mathbf{F}_{\text{even}}$ (respectively, $\mathbf{F}_{\text{odd}}$) stands for the even (respectively, odd) cosine terms in the Fourier series of $\mathbf{F}(a\mathbf{\Psi}(a, \phi) + a\mathbf{B}(a, \phi))$ with respect to $\phi$, Eq. (31) can be split in two equations

$$ 
\mathbf{L}(\Omega^2, \mathbf{\Psi}; \phi) + \frac{1}{a}\mathbf{F}_{\text{odd}}(a, \mathbf{\Psi}, \mathbf{B}; \phi) \\ 
- \mathbf{M}\mathbf{\Psi} \left( \Omega^2 \cos \phi + \frac{1}{2}(\Omega^2)_{,\phi} \sin \phi \right) = 0, \quad (32) $$

$$ \Omega^2 \mathbf{M}\mathbf{B}_{,\phi\phi} + \frac{1}{2}(\Omega^2)_{,\phi} \mathbf{M}\mathbf{B}_{,\phi} + \frac{1}{a} \mathbf{F}_{\text{even}}(a, \mathbf{\Psi}, \mathbf{B}; \phi) = 0, \quad (33) $$

where Eq. (32) (respectively, Eq. (33)) stands for the odd (respectively, even) cosine terms in the Fourier series of Eq. (31).

It should be noticed that $\Omega^2$ and $\mathbf{\Psi}$ are mainly influenced by the odd cosine terms whereas $\mathbf{B}$ is mainly influenced by the even cosine terms. Note also that when $\mathbf{F}$ is an odd function, Eq. (33) is trivially satisfied with $\mathbf{B} = \mathbf{0}$ and Eq. (32) reduces to Eq. (23).

### 4.4. Galerkin procedure to calculate the NNMs

To obtain accurate approximate solutions to Eqs. (32), (33) and (22), a Galerkin method can be implemented. According to the periodic properties, the functions $\Omega^2$, $\mathbf{\Psi}$ and $\mathbf{B}$ can be 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 (34) $$

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

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

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

Inserting Eqs. (34)–(36) into Eqs. (32), (33) and (22), and applying the Galerkin procedure, 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(\Omega_m^2 \cos \phi + \frac{1}{2}(\Omega_m^2)_{,\phi} \sin \phi) \right) \times \cos(2k+1)\phi \, d\phi = 0,
$$

$$ \int_0^{2\pi} \left( \Omega_m^2 \mathbf{M}(\mathbf{B}_m)_{,\phi\phi} + \frac{1}{2}(\Omega_m^2)_{,\phi} \mathbf{M}(\mathbf{B}_m)_{,\phi} + \frac{1}{a} \mathbf{F}_{\text{even}}(a, \mathbf{\Psi}_m, \mathbf{B}_m; \phi) \right) \cos(2k)\phi \, d\phi = 0, $$

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

for $k=0, \dots, m$. This constitutes a set of $2(m+1)n+m+1$ nonlinear equations for the $2(m+1)n+m+1$ unknown coefficients $\Omega_{m,2k}^2$, $\mathbf{\Psi}_{m,2k}$ and $\mathbf{B}_{m,2k}$. The set of algebraic equations can be solved, for given $a$ and $m$ using a Newton–Raphson method, together with an incremental-continuation procedure with respect to the parameters $a$ and/or $m$. The linear normal modes of the underlying linear system can be used as the starting point for $m=0$ and a small value of the amplitude $a$. The order, $m$, of truncated series may be increased to produce more accurate results. As described in Section 3.3, a continuation method (like AMN method) could also be used to solve the set of algebraic equations and hence to detect more easily the bifurcations points.

Once the $\Omega_{m,2k}^2$, $\mathbf{\Psi}_{m,2k}$ and $\mathbf{B}_{m,2k}$ have reached the desired level of accuracy at the desired amplitude level $a$, the differential Eq. (30) can be solved numerically. The resulting time history, $\Phi(t)$, allow us to derive the corresponding time histories of the displacement (29).

### 5. Numerical examples

We consider a 2 d.o.f. nonlinear system composed of a mass $m$ connected to four springs with length $L$ (see Fig. 1). Under the assumption of large displacement, the strain energy can be defined by

$$ W(\mathbf{u}) = \frac{1}{2} \sum_{i=1}^{4} k_i e_i^2 \text{ with } \begin{cases} e_i = \frac{u_i}{L} + \frac{1}{2} \left( \left( \frac{u_1}{L} \right)^2 + \left( \frac{u_2}{L} \right)^2 \right) & \text{for } i=1,2, \\ e_i = - \frac{u_i}{L} + \frac{1}{2} \left( \left( \frac{u_1}{L} \right)^2 + \left( \frac{u_2}{L} \right)^2 \right) & \text{for } i=3,4, \end{cases} $$

where $e_i$ denotes the Green–Lagrange strain of the $i$th spring and $\mathbf{u} = (u_1, u_2)^T$. The equation of motion is then given by

$$ \begin{bmatrix} m & 0 \\ 0 & m \end{bmatrix} \ddot{\mathbf{u}} + \frac{\partial W}{\partial \mathbf{u}}(\mathbf{u}) = 0. \quad (37) $$

Without loss of generality, we will assume that $L=1$ and $m=1$. It should be noted that the Green–Lagrange strain
<figure><img src="image_2.png" /><figcaption>Fig. 1. Two d.o.f. nonlinear system with 4 springs.</figcaption></figure>

$e = \frac{1}{2} \frac{l^2 - L^2}{L^2}$ has been used instead of the more classical engineering strain $e_g = \frac{l-L}{L}$ for the spring. This leads to governing equations with quadratic and cubic nonlinearities as it is the case for beams, plates, and shell. This simple example is in fact representative of thin elastic structure.

## 5.1. Symmetrical case

We consider here the symmetrical configuration of the springs around the mass $m$ defined by $k_1 = k_3$, $k_2 = k_4$, which leads to the following equation of motion:

$$ \ddot{u}_1 + 2k_1u_1 + (k_1 + k_2)u_1(u_1^2 + u_2^2) = 0, \quad (38) $$

$$ \ddot{u}_2 + 2k_2u_2 + (k_1 + k_2)u_2(u_1^2 + u_2^2) = 0. \quad (39) $$

Using the modal representation method it is easy to show that Eqs. (38) and (39) possess two nonlinear normal modes, which can be viewed as the continuation of the normal modes of the underlying linear system. These nonlinear normal modes are defined by the following two pairs of frequency and mode shape functions:

$$ \begin{cases} \Omega_1^2(a, \phi) = 2k_1 + \frac{(k_1+k_2)}{4} a^2 (3 + \cos 2\phi), \\ \Psi_1(a, \phi) = (1, 0)^T, \\ \Omega_2^2(a, \phi) = 2k_2 + \frac{(k_1+k_2)}{4} a^2 (3 + \cos 2\phi), \\ \Psi_2(a, \phi) = (0, 1)^T. \end{cases} $$

Each pair solves the eigenvalue–eigenvector problem (23) and (22). These modes are similar i.e. in the configuration space, the modal line of the first (respectively, the second) mode coincides with the horizontal (respectively, vertical) axis. For each nonlinear normal mode, the motion is given by Eq. (20) solving the differential Eq. (21) with the associated frequency function.

The periodic orbit method has been used with different numbers of time steps. The nonlinear normal modes have been computed by doing these three steps: (i) we have selected a periodic solution with low energy from the linear mode; (ii) we have used a correction algorithm to return on

the nonlinear branch; (iii) we have computed a succession of analytical representations (13) that give a discrete parameterisation in time and continuous parameterisation in energy of the nonlinear normal mode.

For various parameter values $k_1$ and $k_2$, the two methods lead to the same nonlinear normal mode (at least before bifurcation points). Fig. 2 compares the backbone curve in the energy frequency space of the first nonlinear normal mode for $k_1 = 1$ and $k_2 = 2$. The energy is defined by Eq. (2) and the frequency is obtained from the orbit's period. In the modal representation method the backbone

<figure><img src="image_3.png" /><figcaption>Fig. 2. Symmetrical case with $k_1 = 1$ and $k_2 = 2$. Energy versus frequency from the continuation of the first mode: exact solution given by the modal representation method (dashed line), periodic orbit method with $m = 41$ (continuous line) and with $m = 21$ (dotted line).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Symmetrical case with $k_1 = 1$ and $k_2 = 2$. Displacement $u_1$ versus time for $E = 0.734$. Modal representation method (dashed line), periodic orbit method with $m = 41$ (continuous line).</figcaption></figure>
curve has been obtained from Eq. (25), and using the amplitude-energy relationship

$$E(a) = W(a\Psi_1), \qquad (40)$$

deduced from Eq. (2) with the total phase $\Phi$ equal to zero. In this case, the modal representation method gives the exact solution, and, as expected, we can observe that the accuracy of the periodic orbit method increases with the number of time steps used.

Fig. 3 compares periodic motion obtained by the periodic orbit method and the modal representation method for the energy value $E = 0.734$. The periodic solution in the modal representation method has been obtained solving Eq. (21) using the classical Runge-Kutta method with the amplitude value $\bar{a}$ given by $E(\bar{a}) = 0.734$ where $E$ is defined by Eq. (40) and choosing for the initial condition $\phi = 0$ which coincides with the phase condition (6) imposed in the periodic orbit method. In the periodic orbit method the periodic solution has been extracted from the series (13).

## 5.2. Asymmetrical case

We consider now the asymmetrical configuration of the springs around the mass $m$ defined by $k_3 = 0, k_4 = 0$ which leads to the following equations of motion:

$$\ddot{u}_1 + k_1 u_1 + \frac{1}{2} k_1 (3u_1^2 + u_2^2) + k_2 u_1 u_2 + \frac{1}{2} (k_1 + k_2) u_1 (u_1^2 + u_2^2) = 0, \qquad (41)$$

$$\ddot{u}_2 + k_2 u_2 + \frac{1}{2} k_2 (3u_2^2 + u_1^2) + k_1 u_1 u_2 + \frac{1}{2} (k_1 + k_2) u_2 (u_1^2 + u_2^2) = 0 \qquad (42)$$

and we focus here on the first nonlinear mode when $k_1 = 1$ and $k_2 = 2$.

Using the modal representation method, the first nonlinear mode is obtained solving the eigenvalue-eigenvector problem (32), (33) and (22) with the Galerkin procedure described in Section 4.4. The first underlying linear mode ($\omega_1^2 = 1$ and $\psi_1 = (1,0)^T$) has been used as the starting point and $m = 2$ terms have been retained in the finite Fourier series (34)-(36) to approximate the modal functions on the domain $((a, \phi) \in [0,0.36] \times [0,2\pi])$. The two components of the bias vector function $B$ and the two components of the mode shape vector function $\Psi_1$ are shown in Fig. 4. The frequency function $\Omega_1$ is shown in Fig. 5. For $a = 0$, all the modal functions do not depend on the variable $\phi$ and are equal to the corresponding linear mode values. When $a$ increases, the behaviours of the modal functions are monotonous with respect to the variable $a$ and the fluctuations with respect to the variable $\phi$ increase with $a$.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Asymmetrical case with $k_1 = 1$ and $k_2 = 2$. Modal representation of the first nonlinear mode: bias and modal shape vector functions.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Asymmetrical case with $k_1 = 1$ and $k_2 = 2$. Modal representation of the first nonlinear mode: frequency function and invariant manifold in the phase subspace ($U_1, V_1, U_2$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Asymmetrical case with $k_1 = 1$ and $k_2 = 2$. Left: Part of the projection of the invariant manifold in the phase subspace ($u_1, u_2, v_1$) computed by the periodic orbit method. Right: projections of some orbits extracted from the invariant manifold.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Asymmetrical case with $k_1 = 1$ and $k_2 = 2$. Periodic trajectory of the first mode in the phase space for $E = 0.0565$. Modal representation method with $m = 2$ (dashed line), periodic orbit method with $m = 41$ (continuous line).</figcaption></figure>

Thus, the number of terms in the finite Fourier series used in the Galerkin procedure has to be increased to accurately

approximate the modal functions for large values of the amplitude $a$. The softening behaviour of the system is
<figure><img src="image_9.png" /><figcaption>Fig. 8. Asymmetrical case with (a) $k_1 = 1$ and $k_2 = 2$ and (b) $k_1 = 3$ and $k_2 = 1$. Energy versus frequency: from the modal representation method with $m = 2$ (dashed line) and periodic orbit method with $m = 41$ (continuous line).</figcaption></figure>

confirmed in Fig. 5(a). Fig. 5(b) shows a picture of the invariant manifold in the phase subspace ($U_1, V_1, U_2$) where the axis $V_1$ corresponds to the velocity variable $\dot{u}_1$. This three dimensional surface is defined by Eq. (29) and its derivative expression and parameterised by the variable $a$ and $\phi \in [0, 0.36] \times [0, 2\pi]$.

Using the periodic orbit method, the nonlinear normal mode is obtained as in the symmetric case (see previous section). Fig. 6 shows a picture of the invariant manifold of the first nonlinear mode in the phase subspace ($U_1, V_1, U_2$) and some orbits (periodic solutions) extracted from the manifold and projected in the configuration space ($U_1, U_2$) and in the phase subspace ($U_1, V_1$). This approximation of the invariant manifold looks like the approximation obtained by the modal representation method (see Fig. 5). In fact the two approximations are indistinguishable as shown in Fig. 7 where, for example, the trajectory of the periodic orbit for the energy value $E = 0.0565$ obtained by the two methods are plotted in the configuration space ($U_1, U_2$) and in the phase subspace ($U_1, V_1$). With the modal representation method, the periodic orbit has been obtained solving first the scalar differential Eq. (29) with the initial value $\varphi = 0$ and the approximate frequency function $\Omega$ at the amplitude value $\bar{a} = 0.356$ solution of the following algebraic equation:

$$ \frac{\bar{a}^2}{2} (\Psi_{,\phi}(\bar{a}, 0) + \bar{a}\mathbf{B}_{,\phi}(\bar{a}, 0))^T \mathbf{M}(\Psi_{,\phi}(\bar{a}, 0) + \bar{a}\mathbf{B}_{,\phi}(\bar{a}, 0)) + W(\bar{a}\Psi(\bar{a}, 0) + \bar{a}\mathbf{B}(\bar{a}, 0)) = 0.0565, $$

in which the approximated modal functions $\Omega$, $\Psi$ and $\mathbf{B}$ have been used. Finally the mode motion is obtained inserting the trajectory of the total phase $\Phi$ in Eq. (29) evaluated at $a = \bar{a}$.

Fig. 8 shows the backbone curve in the energy frequency space for two pairs of $k_1$ and $k_2$. The two methods have been used, and the results are in good agreement.

## 6. Conclusion

In this paper, two methods have been described for the determination of the nonlinear normal modes of undamped nonlinear mechanical systems. In the periodic orbit method, the nonlinear normal mode is constructed by making the continuation of branches of periodic solutions of the equation of motion. The terms “periodic orbit” means a closed trajectory in the phase space, which is obtained by a suitable energy conserving time integration. In the modal representation approach, the nonlinear normal mode is constructed in terms of amplitude, phase, mode shape, and frequency, with the distinctive feature that the last two quantities are amplitude and total phase dependent. The resulting equation are solved using a Galerkin procedure. Both methods can provide the nonlinear modes at large amplitudes, and they lead to the same results.

These two approaches reveal some differences. The periodic orbit method is more suitable for FEM formulation whereas the modal representation approach is more suitable for parametrical (analytical) description of the modes. Extension to dissipative systems is possible and has been already discussed in [24] but it is an open problem for the periodic orbit method. Finally, the stability analysis of the computed orbit is rather straightforward for the period orbit method using the monodromy matrix [23] whereas it requires additional development for the modal representation approach.

## References

[1] Vakakis AF, Manevitch L, Mikhlin Y, Pilipchuk A, Zevin A. Normal modes and localization in nonlinear systems. New York: John Wiley; 1996.

[2] Gendelman O, Vakakis AF. Transitions from localisation to nonlocalisation in strongly nonlinear damped oscillators. Chaos Solitons Fract 2000;11:1535-42.

[3] Nayfeh AH. Nonlinear interaction: analytical, computational and experimental methods. New York: John Wiley Series in Nonlinear Science; 2000.
[4] Gendelman O, Manevitch LI, Vakakis AF, M'Closkey R. Energy pumping in nonlinear mechanical oscillators: Part I—dynamics of the underlying Hamiltonian systems. ASME J Appl Mech 2001;168: 34-48.

[5] Slater JC, Inman DJ. On the effect of weak non-linearities on linear controllability and observability norms, an invariant manifold approach. ASME J Sound Vib 1997;199(3):417-29.

[6] Pescheck E, Pierre C, Shaw SW. Accurate reduced-order models for a simple rotor blade model using nonlinear normal modes. Math Comput Model 2001;33:1085-97.

[7] Touzé C, Thomas O, Huberdeau A. Asymptotic non-linear normal modes for large-amplitude vibrations of continuous structures. Comput Struct 2004;82:2671-82.

[8] Rosenberg RM. On nonlinear vibrations of systems with many degree of freedom. Adv Appl Mech 1966;242(9):155-242.

[9] Szemplinska-Stupnicka W. The behavior of nonlinear vibrating systems, vols. I and II. Kluwer Academic Publishers; 1990.

[10] Lewandowski R. Computational formulation for periodic vibration of geometrically nonlinear structures, Part 1 and 2. Solids Struct 1996;34(15):1949-64.

[11] Ribeiro P, Petyt M. Nonlinear vibration of plates by the hierarchical finite element and continuation methods. Int J Mech Sci 1999;41: 437-59.

[12] Pérignon F. Vibrations forcées de structures minces, élastiques, non linéaires. PhD thesis, Université de la Méditerranée, Marseille, 2004.

[13] Jézéquel L, Lamarque CH. Analysis of nonlinear dynamical systems by normal form theory. J Sound Vib 1991;149(3):429-48.

[14] Nayfeh AH. Method of normal forms. New York: John Wiley Sons; 1993.

[15] Touzé C, Thomas O, Chaigne A. Hardening/softening behaviour in nonlinear oscillations of structural systems using non-linear normal modes. J Sound Vib 2004;273:77-101.

[16] Nayfeh AH, Mook DT. Nonlinear oscillations. New York: John Wiley & Sons; 1979.

[17] Nayfeh AH, Nayfeh SA. On nonlinear modes of continuous systems. Trans ASME J Vib Acoust 1994;116:129-36.

[18] Shaw SW, Pierre C. Nonlinear normal modes and invariant manifold. J Sound Vib 1991;150:170-5.

[19] Shaw SW, Pierre C. Normal modes for nonlinear vibratory systems. J Sound Vib 1993;164(1):85-124.

[20] Pescheck E, Pierre C, Shaw SW. A new Galerkin-based approach for accurate nonlinear normal modes through invariant manifolds. J Sound Vib 2002;249(5):971-93.

[21] Simo JC, Tarnow N. The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics. Z Angew Math Phys 1992;43:757-92.

[22] Cochelin B, Damil N, Potier-Ferry M. Asymptotic numerical methods and Padé approximants for nonlinear elastic structures. Int J Numer Meth Eng 1994;37:1187-213.

[23] Seydel R. Practical bifurcation and stability analysis, from equilibrium to chaos. 2nd ed. Springer-Verlag; 1994.

[24] Bellizzi S, Bouc R. A new formulation for the existence and calculation of nonlinear mode. In: Bellizzi, Cochelin, Lamarque, editors. Proceedings of EUROMECH 457 (7-9 June 2004, Fréjus, France). Vaulx en Velin: Press of ENTPE; 2004.

[25] Nayfeh AH, Balachandran B. Applied nonlinear dynamics—analytical, computational and experimental methods. New York: John Wiley; 1995.

[26] Van de Vorst EL, Van Campen DH, De Kraker A, Fey RH. Periodic solutions of a multi-dof beam system with impact. J Sound Vib 1996;192(5):913-25.

[27] Doedel EJ. Nonlinear numerics. J Franklin Inst 1997;334B(5/6): 1049-73.

[28] Boutyour EH, Zahrouni H, Poitier-Ferry M, Boudi M. Asymptotic-numerical method for buckling analysis of shell structures with large rotations. J Comput Appl Math 2004;168(1-2):77-85.

[29] Sepulchre JA, MacKay RS. Localised oscillations in conservative or dissipative networks of weakly coupled autonomous oscillators. Nonlinearity 1997;10:679-713.

[30] Galan J, Doedel E, Vanderbauwhede A, Munoz-Almaraz FJ, Freire E. Continuation of periodic orbits in conservative and Hamiltonian systems. Physica D: Nonlinear Phenom 2003;181:1-38.

[31] Crisfield MA. Nonlinear finite element analysis of solids and structures, vol. 1. New York: John Wiley and Sons; 1997.

[32] Zahrouni H, Cochelin B, Potier-Ferry M. Computing finite rotations of shells by an asymptotic numerical method. Comput Meth Appl Mech Eng 1999;175:71-85.

[33] Bellizzi S, Bouc R. A new formulation for the existence and calculation of nonlinear normal modes. J Sound Vib 2005;287(3): 545-69.