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

A new formulation is pres ented here for the exis tence 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 deﬁned by a one-dimens ional nonlinear differential equation governing the total phas e 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 p -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.


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 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 2p-periodic problem, which makes it possible to calculate these quantities using the classical Galerkin procedure in the space of 2p-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.

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 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.
8X; 8Y 2 R n ; ðX À YÞ T ðFðXÞ À FðYÞÞX0: The linear equation where ½qF X ð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 X Ã ðtÞ is a periodic solution with period T40 (i.e.X Ã ðtÞ ¼ X Ã ðt þ TÞ; 8t) of Eq. ( 1), then X ÃÃ ðtÞ ¼ X Ã ðÀtÞ is also a T-periodic solution for Eq. ( 1).Our study will be devoted to the Tperiodic solution of Eq. ( 1) with symmetry X Ã ðtÞ ¼ 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 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 Fð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.

Odd restoring force case
In addition to H1, we assume throughout this section that H2.FðXÞ ¼ ÀFðÀXÞ:

LNMs
We will first deal briefly with the classical linear case in which FðXÞ ¼ ½KX; where ½K is a symmetrical positive definite matrix.The general solution of Eq. ( 1) can be expressed as where ðw l ; o 2 l Þ with o l 40 for l ¼ 1; . . .; n denotes the set of eigenvector-eigenvalue such that ½Kw l ¼ ½Mw l o 2 l ; w T l ½Mw l ¼ 1 and w T l ½Mw k ¼ 0; lak: The pairs ða l ; j l Þ 2 R þ Â ½0; 2p of constant real numbers set the initial conditions.It follows that Eq. ( 1) possesses n-periodic solutions with period T l ¼ 2p=o l ; l ¼ 1; . . .; n; respectively, which can be written as (index l is omitted) XðtÞ ¼ wa cos FðtÞ; (5) The scalar function Fð:Þ 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.

NNMs
As suggested by the linear case, a synchronous periodic solution to Eq. ( 1) is sought in the form XðtÞ ¼ Wða; FðtÞÞa cos FðtÞ;

_
FðtÞ ¼ Oða; FðtÞÞ; Fð0Þ ¼ j: 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 Fð:Þ: Four quantities are involved in Eqs.(7) and (8): two given scalar quantities, the amplitude ðaX0Þ and the phase (j 2 ½0; 2p) which set the initial conditions of the vibration in NNM motion on the configuration space.Two unknown functions, W and O; 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 W is required.In this study, we will adopt (without lost of generality) the following condition: The modal vector W and the frequency O will be viewed here as amplitude and total phase dependent, and denoted W ¼ Wða; FÞ; O ¼ Oða; FÞ; and will be searched for as 2p-periodic functions with respect to F for fixed a.More specifically, according to the symmetry condition H2, these functions will be searched for as even periodic functions with respect to F with period 2p; that is with Fourier series Oða; FÞ ¼ X define a ''synchronous'' periodic oscillation [2], where all the material points in the system pass through zero ( simultaneously.The set, in the phase space, defined by Eqs. ( 7) and (11) taking the initial conditions a and j 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] XðtÞ ¼ aWðaÞ cosðOðaÞt þ jÞ; 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.

A periodic eigenvalue-eigenfunction problem
The objective now is to characterize the pair ðO; WÞ which defines an NNM.Substituting Eq. (7) into Eq.(1), using Eq. ( 8 where F looks like an independent variable.This variable will therefore be denoted f in what follows. Differentiating with respect to f; using the differential rule 2OqO=qf ¼ qO 2 =qf; Eq. ( 13) becomes where the differential operator L is given by LðO The frequency function, O; appears in Eq. ( 14) by its square, O 2 : So, in the following, we will look for the pair ðO 2 ; WÞ: we must solve, for the fixed parameter a, the 2p-periodic nonlinear eigenvalue ðO 2 Þ-eigenvector ðWÞ problem defined by Eqs. ( 14) and ( 16).These equations reduce to ½KW ¼ ½MWO 2 in the linear case.We will focus on periodic solutions ðO 2 ða; :Þ; Wða; :ÞÞ which satisfy the following properties: Moreover, if we assume that the eigenvalues associated with the pair of matrices ð½M; ½q 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 ðo 2 l ; w l Þ where ðo 2 l ; w l Þ denotes the corresponding normal mode of the underlying linear system (3).
Furthermore, as a !0; the following limits hold: Consequently there exist n NNMs which can be viewed as a continuation of the LNMs of the underlying linear system (3).

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 where X i denotes the ith component of X: This leads to the following statement: an NNM is similar if and only if the ½M-normalized vector W does not depend on F: In this case, Eq. ( 14), using Eq. ( 16), reduces to and substituting Eq. (20) into Eq.( 14) yields FðaW cos fÞ ¼ ðW T FðaW cos fÞÞ½MW: This leads to the following statement: the nonlinear system (1) admits an SNNM if and only if Eq. ( 21) possesses a ½M-normalized f-independent vector solution W: The associated frequency is given by This expression defines a regular, even, 2p-periodic function with respect to f with a Fourier series of the form O 2 ða; fÞ ¼ P 1 k¼0 O 2 2k ðaÞ cos 2kf: Comment: It can be easily checked that any function of the following form yields SNNMs: where ½K and ½K 1 are symmetrical and definite positive matrices such that the eigenvalue problem associated with the pair of matrices ð½M; ½KÞ and ð½M; ½K 1 Þ leads to the same eigenvectors, and f ð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); FðXÞ ¼ a½MX þ GðXÞ where aX0 and G i ðXÞ; i ¼ 1; 2; . . .; n; are homogeneous polynomials of degree 2k þ 1: where ½K and ½C are real square matrices, ½K is symmetrical and positive definite and l is a real number.Eq. ( 21) reduces to Hence the mode shapes of the SNNMs coincide with those of the underlying linear system ðl ¼ 0Þ: Let ðo p ; w p Þ be a normal mode of the linear system.Using w T p ½Mw p ¼ 1 and w T p ½Kw p ¼ o 2 p ; Eq. ( 22) gives The pair ðO 2 p ; w p Þ completely defines an SNNM which tends towards ðo 2 p ; w p Þ as a tends towards zero.The normal mode can be a softening or hardening mode depending on the sign of lðw T p ½Cw p Þ (provided hypothesis H1 is fulfilled).The corresponding period of the oscillations is given by ( 10) where Kð:Þ denotes the elliptic function of the first kind.Note that Tða ¼ 0Þ ¼ 2p=o 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.

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 ð2pÞ: To do this, a Galerkin method is implemented [10].The unknown functions (square-frequency O 2 and mode shape C), are expanded into a finite Fourier series with respect to the variable f according to W m;2k ðaÞ cos 2kf: The determining equations read, for k ¼ 0; . . .
m;2k and W 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 f; the following equations 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 incrementalcontinuation procedure with respect to the parameters a and/or m.For a ¼ a þ Da 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; W T 0;0 ½MW 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 and/or this sup-norm.
Once the O 2 m;2k 's and W 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, FðtÞ; allow us to derive the corresponding time histories of the displacement (7) and velocity (11).

Example 1
The first example is taken from Ref.
where o 1 ¼ 0:689 and o 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 O 2 0;0 ¼ 0:689 2 ; W 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 Wða; fÞ and the frequency Oða; fÞ 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 Wða; 0Þ go beyond the range of small motions, and large-amplitude motion is generated.At fixed f; the frequency O increases with the amplitude a.The maximum values of O in the cross-sections defined by a ¼ C t are obtained at f ¼ 0 and p 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 _ 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 j ¼ 0 and Eq. ( 7) is performed using the frequency Oð1:88; FðtÞÞ and mode shape Wð1:88; FðtÞÞ with (m ¼ 6)-harmonics, while on the other hand, a direct numerical integration of Eq. ( 1) with the corresponding initial conditions ð1:88C 1 ð1:88; 0Þ; 0; 1:88C 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 f (Fig. 4(a)).
Similar results were obtained for the second NNM as a continuation of the corresponding linearized mode O 2 0;0 ¼ 3:244 2 ; W 0;0 ¼ ð0; 1Þ: The m ¼ 3 order approximation was sufficiently accurate to describe its behaviour.

General restoring force case
Throughout this section, we only assume H1.

NNMs
In the case where FðXÞ is not necessarily an odd function, the NNMs take the following form: where B denotes a vector function.
Five quantities are involved in Eqs. ( 30) and ( 31): two given scalar quantities, the amplitude ða40Þ and the phase (j 2 ½0; 2p), which set the initial conditions of the vibration in NNM motion on the configuration space.Three unknown functions, W; O and 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 F with period p; that is with Fourier series 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 has been considered, ignoring the phase dependence on the frequency, mode shape and bias term.

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 where m denotes the order of the truncated series.
The corresponding determining equations read ; and O 2 m;2k : As previously, an incremental-continuation numerical procedure based on the Newton-Raphson algorithm can be implemented.

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 ; where denotes the strain-displacement of the rth spring, the equations of motion can be written, assuming L ¼ 1; as follows: where 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: o 1 ¼ 1 and o 2 ¼ 2; reflecting the presence of a two-order internal resonance, o 2 ¼ 2o 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; FÞ 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][14][15][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 f ¼ 0) and the direct numerical integration of Eqs. ( 41) and (42) are indistinguishable.

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 2p-periodic differential problem.As in the linear case, the eigenvalue (frequency) is always positive and can be expressed in terms of the eigenvector (mode shape).The eigenvectoreigenvalue 2p-periodic problem is numerically solved using the classical Galerkin procedure in the space of 2p-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 where the pair ðW l ; O l Þ stands for the l-(coupled) mode, and a ¼ ða 1 ; a 2 ; . . .; a n Þ T ; U ¼ ðF 1 ; F 2 ; . . .; F 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.