# Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model

Aurélien Grolet, Fabrice Thouverez

► To cite this version:

Aurélien Grolet, Fabrice Thouverez. Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model. Journal of Sound and Vibration, 2012, 331 (12), pp.2911-2928. 10.1016/j.jsv.2012.02.008. hal-02424287

HAL Id: hal-02424287

https://hal.science/hal-02424287

Submitted on 27 Dec 2019

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
Free and forced vibration analysis of a nonlinear system
with cyclic symmetry: Application to a simplified model

Aurelien Grolet*, Fabrice Thouverez

Ecole Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes, 36, avenue Guy de Collongue, 69134 Ecully Cedex, France

This work program is devoted to studying the nonlinear dynamics of a structure with cyclic symmetry under conditions of geometric nonlinearity, through the use of the harmonic balance method (HBM). In order to study the influence of nonlinearity due to the large deflection of blades, a simplified model has been developed. This approach leads to a system of linearly coupled, second-order nonlinear differential equations, in which nonlinearity appears via cubic terms. Periodic solutions, in both the free and forced cases, are sought by applying HBM coupled with an arc-length continuation method. Solution stability has been investigated using Floquet's theorem. In addition to featuring similar and nonsimilar nonlinear modes, the unforced system is known to contain localized nonlinear modes that arise from branching point bifurcation at certain vibration amplitudes. In the forced case, these nonlinear modes give rise to a complex dynamic behavior. Many bifurcations can take place, thus leading to strong or weak localization that may or may not be stable. In this study, special attention has been paid to the influence of excitation on dynamic responses. Several cases of excitation have been analyzed herein: localized excitation, and low-engine-order excitation. In the case of low-engine-order excitation, sensitivity of the response to a perturbation of this excitation type has been investigated, and it has been shown that for a localized, or sufficiently detuned excitation, several solutions can coexist, some of which are represented by closed curves in the Frequency-Amplitude domain. These various solutions overlap when increasing the force amplitude, leading to forced nonlinear localization. Because closed curves are not tied up with the basic nonlinear solution, they can easily be overlooked. In this study, they have been calculated using a sequential continuation with the force amplitude as a parameter.

# 1. Introduction

This paper is intended to study both free and forced nonlinear vibrations of a structure with cyclic symmetry [1-3], composed of *n* identical substructures subjected to large displacements. This type of system appears, in particular, in the context of modeling bladed wheels and spatial antennas [3,4]. The nonlinearities involved stem from the large displacements of the substructures (i.e. geometric type nonlinearities). Establishing the equations for such structures leads to a system of second-order, coupled differential equations containing polynomial-type nonlinear terms.

In the linear case, due to perfect problem symmetry, a study of these equations reveals a majority of double eigenfrequencies, corresponding to distinct mode shapes. These mode shapes are associated with diameter vibration modes; for a given eigenfrequency, the associated mode shapes feature mutually orthogonal diameters. Moreover, for weakly coupled

<footnote>* Principal corresponding author.</footnote>

<footnote>E-mail addresses: aurelien.grolet@ec-lyon.fr, aurelien.grolet@ec2009.ec-lyon.fr (A. Grolet), fabrice.thouverez@ec-lyon.fr (F. Thouverez).</footnote>
and slightly detuned (linear) systems (i.e. systems with minor imperfections that lead to symmetry breaks), it has been demonstrated that localization phenomena may arise [5]. Such phenomena produce confined motion on one or more system substructures, leaving the others practically motionless.

In the nonlinear case, Vakakis [6,4,7] showed that localization phenomena are indeed possible, in either a free or forced mode, for a weak enough coupling between substructures or a strong enough nonlinearity. This author also studied the transition between localization and non-localization as a function of the coupling/nonlinearities ratio; he found that localization was only possible for small values of this ratio and in cases where bifurcations appear as this ratio increases, thus leading to a disappearance of localized solutions. Localization phenomena have also been identified by many studies such as [8–10].

During more recent research [11,3], Peeters applied a shooting method to calculate the nonlinear normal modes of a system with cyclic symmetry, which has the effect of exposing both the similar and nonsimilar normal modes and even localization phenomena for some nonlinear normal modes. He showed that mode localization arise from bifurcations of the free undamped systems.

In the scope of the present paper, we propose studying both the free and forced responses of a simplified model of a bladed wheel. Contrary to Peeters, we will use the harmonic balance method (HBM) to derive solutions. This technique will yield, as part of its very nature, analytical as well as numerical studies. We aim at conducting bifurcation analysis (whether in the free or forced case) of the obtained HBM equations and thus confirm that mode localization arise from bifurcation of the free undamped system. Moreover, in the forced case, we will focus on localization effects and the links between disymmetric loading and localized responses. Special attention will be paid to the influence exerted by loading level on the forced response, which allows drawing original observations on the appearance of periodic solutions displayed in the form of closed curves in the frequency-amplitude diagram. Stability analysis will also be conducted, which will reveal that some parts of the closed curve solution may be stable, thus implying that this new phenomenology may be important in the design of bladed disks.

## 2. Problem presentation and equations of motion

Let us consider a structure with cyclic symmetry composed of *n* identical substructures. It is decided to model these sectors by means of beams embedded onto a fixed foundation at one end and axially restricted on the other end. The coupling between beams is introduced via a linear stiffness running between two consecutive substructure (Fig. 1). When undergoing relatively large displacements, substructure can be subjected to geometric nonlinearities. Many studies have been conducted to examine the effect of geometric nonlinearities on thin structures (e.g. [12,13]).

Given the level of problem symmetry, the initial focus will involve just one beam; afterward, the other beams will be integrated along with the linear coupling stiffening elements.

### 2.1. Derivation of strain and kinetic energies for a single beam

Now assume a beam of length *L*, constant section *S* and constant quadratic momentum *l*, associated with a Cartesian coordinate system (*O*x, *y*,*z*), embedded at the end (*x* = 0) and axially restricted at the end *x*=$L$ (Fig. 2). The displacement of a single coordinate point (*x*,*y*) in direction (*O*x) (resp. (*O*y*)) is denoted $\tilde{u}$ (respectively $\tilde{v}$).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Cyclic system used for establishing equations.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Diagram of the beam considered and boundary conditions.</figcaption></figure>
Boundary conditions for such a beam are given by the following:

* Kinematic boundary conditions:

$$ \tilde{u}(0,y) = \tilde{u}(L,y) = \tilde{v}(0) = \frac{\partial \tilde{v}}{\partial x}(0) = 0 \quad (1) $$

* Natural boundary conditions:

$$ M(L) = T(L) = 0 \quad (2) $$

where M is the bending moment and T the shear force.

In the following we will derive an expression for the different energies (strain, kinetic and work of external forces) of a single beam in terms of the transverse and axial displacements. Then we will show that applying Hamilton's principle and neglecting axial inertia allow to remove the axial displacement from the different energies thus leading to energies expressed only in terms of the transverse displacement. After introducing cyclic symmetry, these energies will be used to derive a discretization of the motion's equations through a Rayleigh-Ritz approach along with Lagrange's equations.

**Strain energy:** In the following discussion the Lagrangian formulation is adopted along with classical Von Karman theory hypothesis [14]. Among them, the Kirchhoff hypothesis is introduced so that stresses in the direction normal to the beam middle line are negligible, and strains vary linearly within the beam thickness. In this case, displacements of the assembly are entirely determined by middle line displacements u and w, by the following relation:

$$ \mathbf{u}(x,y) = \begin{pmatrix} \tilde{u}(x,y) \\ \tilde{v}(x) \end{pmatrix} = \begin{pmatrix} u(x)-y\frac{\partial u}{\partial x} \\ v(x) \end{pmatrix} \quad (3) $$

The nonlinear relation between displacements u, v and strains of the middle line $\epsilon$ is given by the Von Karman relationship for beams i.e.:

$$ \epsilon = \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial v}{\partial x} \right)^2 \quad (4) $$

The total strain $\tilde{\epsilon}$ is expressed as

$$ \tilde{\epsilon} = \epsilon - yk_z \quad (5) $$

where $k_z = \partial^2 v / \partial x^2$ is the change of curvature.

The elastic strain energy of the beam is given by

$$ U = \frac{1}{2} \int_{0}^{L} \left( \int_{\Sigma} \sigma \tilde{\epsilon} d\Sigma \right) dx \quad (6) $$

where $\sigma$ is the Kirchhoff stresses, related to Green's strain $\tilde{\epsilon}$ by a standard Hooke's law given by

$$ \sigma = E\tilde{\epsilon} \quad (7) $$

where E is Young's modulus of the material.

Finally by substituting Eq. (7) into Eq. (6) and doing integration, we obtain the following expression for the strain energy:

$$ U = \frac{1}{2} \int_{0}^{L} \left( \int_{\Sigma} E \epsilon^2 d\Sigma \right) dx + \frac{1}{2} E I \int_{0}^{L} \left( \frac{\partial^2 v}{\partial x^2} \right)^2 dx \quad (8) $$

**Kinetic energy:** If the inertia of rotation is neglected the kinetic energy of a beam is given as

$$ T = \frac{1}{2} \rho \int_{0}^{L} \left( \int_{\Sigma} (\dot{u}^2 + \dot{v}^2) d\Sigma \right) dx \quad (9) $$

where $\rho$ is the density of the material.

**External force:** In the following discussion, only the case of a punctual periodic force orthogonal to the beam applied in $x_f = L$ will be considered. The work of such a force is given by the following relation:

$$ W = v(L)g(t) \quad (10) $$

where $g(t)$ is a periodic function.

**Hamilton's principle:** The previous developments (Eqs. (8)-(10)) allow to derive an expression for the Lagrangian of a single beam system. The application of Hamilton's principle gives the following continuous equation of motion for the
axial displacement u:

$$ \rho S \ddot{u} + \frac{\partial}{\partial x} \left[ ES \left( \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial v}{\partial x} \right)^2 \right) \right] = 0 \quad (11) $$

By defining the tension along the beam as $N = \int_{\Sigma} \sigma d\Sigma = ES\epsilon$ and by neglecting axial inertia effects, Eq. (11) can be rewritten as

$$ \frac{\partial N}{\partial x} = 0 \quad (12) $$

From Eq. (12) we can conclude that $N$ is constant along the beam. By integrating $N$ from 0 to $L$ we have

$$ \int_{0}^{L} N \, dx = LN = ES \int_{0}^{L} \left( \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial v}{\partial x} \right)^2 \right) dx \quad (13) $$

Considering the kinematic boundary conditions in Eq. (1) (i.e. $u(0) = u(L) = 0$) gives $\int_{0}^{L} (\partial u / \partial x) dx = 0$ and the in plane displacement $u$ vanishes from the expression of the normal stress, which is ultimately given by the following relation:

$$ N = \frac{ES}{L} \int_{0}^{L} \left( \frac{\partial v}{\partial x} \right)^2 dx \quad (14) $$

We then deduce the following expression for the in plane strain $\epsilon$:

$$ \epsilon = \frac{N}{ES} = \frac{1}{2L} \int_{0}^{L} \left( \frac{\partial v}{\partial x} \right)^2 dx \quad (15) $$

Substituting Eq. (15) into Eq. (8) results in an expression for the strain energy $U$ only in terms of the transverse displacement $v$ given by

$$ U = \frac{ES}{8L^2} \int_{0}^{L} \left( \int_{0}^{L} \left( \frac{\partial v}{\partial x} \right)^2 dx \right)^2 dX + \frac{1}{2} EI \int_{0}^{L} \left( \frac{\partial^2 v}{\partial x^2} \right)^2 dx \quad (16) $$

Finally, by neglecting the axial inertia the kinetic energy in Eq. (9) is reduced to the following:

$$ T = \frac{1}{2} \rho \int_{0}^{L} \left( \int_{\Sigma} \dot{v}^2 d\Sigma \right) dx \quad (17) $$

In the remaining of the paper the strain and kinetic energies will be defined by Eqs. (16) and (17) and will be used for the computation of the Lagrangian of the system. But before that we will present how we switch from a single beam to the cyclic symmetry model presented in Fig. 1.

## 2.2. Incorporation of cyclic symmetry

Let us return to the case of the structure with cyclic symmetry composed of $n$ identical sectors. For $1 \le j \le n$, $v^j$ denotes the transverse displacement of the $j$th beam relative to its coordinate system, $U^j$ its strain energy, $T^j$ its kinetic energy and $W^j$ the work of forces being applied on this beam. The expressions for $U^j$, $T^j$ and $W^j$ have been obtained by replacing $v$ by $v^j$ in Eqs. (16), (17), (10). The expressions of strain energy $U$, kinetic energy $T$, and work $W$ for the complete system can then be given by

$$ U = \sum_{j=1}^{n} U^j, \quad T = \sum_{j=1}^{n} T^j \quad \text{and} \quad W = \sum_{j=1}^{n} W^j \quad (18) $$

Let us now assume that the coupling between beams is modeled by a set of massless linear springs, all with an identical constant of stiffness $k$. The springs are set at points in $x_r = L/4$ on all beams. The energy stored by such a spring located between beams $j$ and $j+1$ is expressed by

$$ V^j = \frac{1}{2} k \left( v^j \left( \frac{L}{4} \right) - v^{j+1} \left( \frac{L}{4} \right) \right)^2 \quad (19) $$

for $1 \le j \le n$ according to the convention $j+1 = 1$ if $j = n$ (cyclic symmetry). The total energy $V$ stored by the $n$ springs can then be written as

$$ V = \sum_{j=1}^{n} V^j \quad (20) $$
## 2.3. Discretization and equations of motion

Now, let us introduce the discretization protocol for transverse displacement $v^j$. This discretization routine is performed by relying on a Rayleigh-Ritz approach [15], which generates an approximated model of the structure. The transverse displacements ($v^j$) are interpolated using the following expression:

$$v^j(x,y,t) = \sum_{i=0}^{N} \lambda_i^j(t) \Phi_i(x,y) \qquad (21)$$

where the ($\Phi_i$) are kinematically admissible functions, $\lambda_i^j$ the participation of these functions in the displacement $w^j$, and $N$ the number of Ritz functions selected for the interpolation step. By using this method, the transition is made from a continuous problem to a discrete one of dimension $nN$ where the unknowns are the variables ($\lambda_i^j$)$_{1 \le i \le N, 1 \le j \le n}$. By substituting $v^j$ into Eqs. (19) and (21), an expression of the various energies can be derived as a function of ($\lambda_i^j$).

The motion equations are then provided by the Lagrange equations, written as follows:

$$\frac{d}{dt} \frac{\partial T}{\partial \dot{\lambda}_i^j} + \frac{\partial U}{\partial \lambda_i^j} + \frac{\partial V}{\partial \dot{\lambda}_i^j} = \frac{\partial W}{\partial \lambda_i^j} \quad \text{for } 1 \le i \le N, 1 \le j \le n \qquad (22)$$

For numerical purpose, we arbitrarily choose to set $n=6$ (only six substructures) and $N=1$ (only one Ritz function) thus yielding a nonlinear problem with 6 degrees of freedom that will serve as an example for the present work. These numerical values are sufficient to illustrate the main effects, but the procedure remains the same for any values of $n$ and $N$.

The single Ritz function selected for this interpolation step is $\Phi = (x/L_x)^2$, which verifies the kinematic boundary conditions given in Eq. (1). Applying Lagrange's equations (Eq. (22)) for the discretization $v^j = \lambda_j(x/L)^2$ results in the following equation of motion for the jth beam:

$$\ddot{\lambda}_j + (\alpha + 2c)\dot{\lambda}_j + \beta\dot{\lambda}_j^3 - c\lambda_{j+1} - c\lambda_{j-1} = f_j(t) \qquad (23)$$

with

$$\alpha = \frac{20EI}{\rho SL^4}, \quad \beta = \frac{40ES}{9\rho SL^4}, \quad c = \frac{5k}{256\rho SL}, \quad f_j = \frac{5g_j}{\rho SL} \qquad (24)$$

By setting vector $X = (\lambda_j)_{1 \le j \le 6}$ and the convention $X^3 = (\lambda_j^3)_{1 \le j \le 6}$, the matrix form of Eq. (23) is given by the following:

$$[\mathbf{M}]\ddot{\mathbf{X}} + [\mathbf{K}]\mathbf{X} + \beta \mathbf{X}^3 = \mathbf{F}(t) \qquad (25)$$

with $[\mathbf{M}] = [\mathbf{I}]_6$ and

$$[\mathbf{K}] = \begin{pmatrix} \alpha + 2c & -c & 0 & 0 & 0 & -c \\ -c & \alpha + 2c & -c & 0 & 0 & 0 \\ 0 & -c & \alpha + 2c & -c & 0 & 0 \\ 0 & 0 & -c & \alpha + 2c & -c & 0 \\ 0 & 0 & 0 & -c & \alpha + 2c & -c \\ -c & 0 & 0 & 0 & -c & \alpha + 2c \end{pmatrix} \qquad (26)$$

The numerical values to be input throughout the remainder of this study are given by the following:

$$L = 1.5 \text{ m}$$

$$\rho SL = 105.3 \text{ kg}$$

$$\frac{ES}{L} = 2.484 \times 10^9 \text{ N m}^{-1}$$

$$\frac{EI}{L^3} = 4.6123 \times 10^4 \text{ N m}^{-1}$$

$$k = 8 \times 10^5 \text{ N m}^{-1} \qquad (27)$$

which in turn correspond to the parameter values in Eq. (23):

$$\alpha = 8,7662 \times 10^3 \text{ s}^{-2}, \quad c = 148.36 \text{ s}^{-2}, \quad \beta = 4,6752 \times 10^7 \text{ m}^{-2} \text{ s}^{-2} \qquad (28)$$
### 3. Solution method: harmonic balance method (HBM)

#### 3.1. Principle of the HBM method

The harmonic balance method is a widely used technique for studying nonlinear systems, with a large number of application illustrations available. This method offers the ability to treat the case of a highly nonlinear system with, for example, friction [16,17] or geometric nonlinearities [18,19]. The following system of nonlinear differential equations will now be considered:

$$[\mathbf{M}]\ddot{\mathbf{X}} + [\mathbf{K}]\mathbf{X} + \mathbf{K}_{nl}(\mathbf{X}) = \mathbf{F}(t) \qquad (29)$$

where $\mathbf{X}$ is a vector of unknowns with dimension $n$, $[\mathbf{M}]$ and $[\mathbf{K}]$ are the system's mass and stiffness matrices respectively; $\mathbf{K}_{nl}$ is the nonlinear forces and $\mathbf{F}$ the external loading which is assumed to be periodic with period $T = 2\pi/\omega$.

The harmonic balance method consists of identifying periodic solutions $\mathbf{X}$ to the system in Eq. (29) in the form of a truncated Fourier series, i.e.:

$$\mathbf{X}(t) = \mathbf{A}_0 + \sum_{k=1}^{N_h} \mathbf{A}_k \cos(k\omega t) + \mathbf{B}_k \sin(k\omega t) \qquad (30)$$

By including this development in Eq. (29), and by projecting equations on the Fourier basis, as follows:

$$\int_{0}^{2\pi/\omega} ([\mathbf{M}]\ddot{\mathbf{X}} + [\mathbf{K}]\mathbf{X} + \mathbf{K}_{nl}(\mathbf{X})) \cos(k\omega t) dt = \int_{0}^{2\pi/\omega} \mathbf{F}(t) \cos(k\omega t) dt \qquad (31)$$

$$\int_{0}^{2\pi/\omega} ([\mathbf{M}]\ddot{\mathbf{X}} + [\mathbf{K}]\mathbf{X} + \mathbf{K}_{nl}(\mathbf{X})) \sin(k\omega t) dt = \int_{0}^{2\pi/\omega} \mathbf{F}(t) \sin(k\omega t) dt \qquad (32)$$

then a system of $n(2N_h+1)$ nonlinear algebraic equations is obtained, containing $n(2N_h+1)$ unknowns $(\mathbf{A}_k)_{0 \le k \le N_h}$, $(\mathbf{B}_k)_{1 \le k \le N_h}$. Note that in the case of a free undamped system, the angular frequency $\omega$ is not known a priori, and we have to consider it as a variable.

One inherent advantage of this method is its absence of a hypothesis regarding the order of magnitude of nonlinearities: the same treatment procedure is employed regardless of whether a given system is strongly or weakly nonlinear. The most important parameter of this method is the number of selected harmonics $N_h$; this number in fact is not typically known ahead of time, a situation that necessitates convergence studies in order to ensure accurate solution estimation. On the other hand, in the case where the number of harmonics selected is high, the solution procedure can quickly become complex and difficult. In many cases however, it turns out that the solution rapidly converges in terms of harmonics, which leads to systems with more practical dimensions.

#### 3.2. Numerical procedure

Generally speaking, the system resulting from the application of the HBM can be rewritten in the following form:

$$\mathbf{G}(\mathbf{Y}, \omega) = 0 \qquad (33)$$

where $\mathbf{Y}$ is a vector of dimension $n(2N_h+1)$ that corresponds to the HBM unknowns $(\mathbf{A}_k)_{0 \le k \le N_h}$ and $(\mathbf{B}_k)_{1 \le k \le N_h}$, and $\mathbf{G}$ a nonlinear value function within a space of dimension $n(2N_h+1)$.

The solution to Eq. (33) will make use of a continuation technique [11,20–22] that consists of numerically searching for a series of $M+1$ points $(\mathbf{Y}^i, \omega^i)_{0 \le i \le M}$ that verifies the following:

$$\mathbf{G}(\mathbf{Y}^i, \omega^i) \le \epsilon \qquad (34)$$

where $\epsilon$ is a small parameter introduced to determine the level of solution accuracy.

The continuation procedure takes place in two stages: a prediction phase, followed by a correction phase. The predictor yields an estimation of the solution sought based on the previously determined solution points, while the corrector iteratively performs corrections on the predicted solution so as to derive a solution that satisfies criterion in Eq. (34).

An arc length type of continuation will be used herein [21]; as such, the solutions are to be parameterized by the curvilinear abscissa $s$, with $\omega$ being considered as a variable and with the solutions to Eq. (33) being of the form $(\mathbf{Y}(s), \omega(s))$.

#### 3.3. Solution stability and bifurcations

To determine the stability of periodic solutions, we turn to the well known Floquet's theory, and we compute the so-called monodromy matrix. The eigenvalues of the monodromy matrix offer the ability to draw a conclusion on both the stability of solutions and the nature of bifurcations likely to appear [20].

The focus now turns to two types of bifurcation (of the nonlinear algebraic system) that very often appear in this study: turning point type bifurcation, and branching point type bifurcation. For the turning point bifurcation, no special procedure is applied since the arc length continuation algorithm is adapted to handle such situations.
The detection of bifurcation points is done by monitoring the determinant $|J|_{[y]}$ of the Jacobian matrix $J|_{[y]} = \partial G / \partial Y$ of the system in Eq. (33): when this determinant is zero at some point, then the matrix $J|_{[y]}$ is singular and this point is a bifurcation point. In order to precisely obtain the bifurcation point, we used an algorithm described in [20] that finds the solution $(Y_b, \omega_b)$ such that

$$G(Y_b, \omega_b) = 0$$

$$\frac{\partial G}{\partial Y}(Y_b, \omega_b)T = 0$$

$$\|T\| = 1 \qquad (35)$$

Once bifurcation points have been obtained, the type of bifurcation is determined by the rank of the $M_h \times (M_h+1)$ matrix $J = [J|_{[y]} \ J|_{[\omega]}]$ (where $J|_{[\omega]} = \partial G / \partial \omega$). At the level of the branching points, the matrix $J$ has a rank of at most $M_h-1$. Several linearly independent tangents $T_i$ can then be found such that $J|_{[y]}T_i = 0$, which in practice corresponds to several prediction directions for the continuation algorithm. In this study, we computed the prediction direction by calculating the eigenvectors $\Psi_i$ of the matrix $J|_{[y]}$ associated with the zero eigenvalue at the bifurcation point [23].

Another way of computing bifurcated branches is the so-called perturbation method which consists of adding a small constant perturbation to the function $G(Y, \omega)$ (Eq. (33)) in the vicinity of the bifurcation point. This step serves to separate the structurally unstable branching points. A continuation can thus be carried out on a few points with the perturbed system, during the time it takes to extend beyond the branching point and lock in on another solution [20].

## 4. Free response of a 6 degree-of-freedom system: nonlinear normal modes

The objective in this section is to estimate the free periodic solutions of the system in Eq. (25). The natural approach here would thus be to make use of nonlinear normal modes (NNM) which have already yielded an abundance of scientific output. In [24,25], Vakakis and Kerschen conducted a literature review of results relative to these NNM, which may be defined in one of the several ways. The most widespread definition has been provided by Rosenberg [26] who defined NNM as vibrations in unison of the free undamped system:

$$[M]\ddot{X} + [K]X + \beta X^3 = 0 \qquad (36)$$

This definition imposes that the degrees of freedoms vibrate with the same frequency, reaching their extreme values while simultaneously passing through zero. Consequently, all displacements may be expressed as a function of a single (reference) degree of freedom. This parameterization allows calculating and representing the NNM by curves within the configuration space. When these curves are straight lines (i.e. a linear relation between the degree of freedom (dof) and the reference coordinate), the NNM are said to be similar and their mode shapes are independent on vibration amplitude. In the general case, the modal lines are not straight-line curves and the corresponding NNM are said to be nonsimilar, since both their mode shapes and eigenfrequencies vary with respect to vibration amplitude. These modal curves are tangent to the straight lines that correspond to the eigenmode of the linear system for small vibration amplitudes. Among the other approaches available for defining NNM, let us highlight the method of invariant manifolds proposed by Shaw and Pierre [27], the multiple scales method [20,4], the normal form method [28], the harmonic balance method [19], and the shooting method [11].

We will start by studying the linear system underlying the problem expressed in Eq. (25) which will serve as a starting point for the nonlinear analysis. The next step will concentrate on identifying system NNM through the use of the harmonic balance method described in Section 3. An analytical assessment will first be drawn before numerically searching for solutions by means of coupling the HBM with an arc-length continuation.

### 4.1. Study of the underlying linear system

Let us consider the free undamped linear system underlying the problem displayed in Eq. (25):

$$[M]\ddot{X} + [K]X = 0 \qquad (37)$$

The linear eigenmodes (both eigenfrequencies and mode shapes) of Eq. (37) are given by the solution to the following eigenvalue equation:

$$[K]\Phi = \omega^2[M]\Phi \qquad (38)$$

Matrix $[K]$ contains simple as well as double eigenvalues (we recall that in this study $[M]$ equals the identity matrix), yet the double eigenvalues are predominant (i.e. a simple eigenvalue if the dimension of $[K]$ is odd-numbered, versus two simple eigenvalues if the dimension is even). Each double eigenvalue corresponds to two distinct eigenvectors that create a two-dimensional solution space. In the general case, the eigenvalues of matrix $[K]$ are given by the formula derived from [1]:

$$\omega_p^2 = \alpha + 2c(1 - \cos(\theta_p)) \quad \text{with} \quad \theta_p = \frac{2p\pi}{n} \qquad (39)$$
with $p$ being the mode diameter number, such that $0 \le p \le (n-1)/2$ if $n$ is odd and $0 \le p \le n/2$ if $n$ is even. For a weak coupling coefficient $c$, the eigenfrequencies lie relatively close to one another. The linear mode shapes are given by

$$ \Phi_0 = [1, 1, 1, 1, 1, 1]^T $$

$$ \Phi_p^c = [\cos(\theta_p), \cos(2\theta_p), \dots, \cos(n\theta_p)]^T $$

$$ \Phi_p^s = [\sin(\theta_p), \sin(2\theta_p), \dots, \sin(n\theta_p)]^T \quad (40) $$

In our specific case, application of these formulae yields the eigenfrequencies (expressed in rad s⁻¹) along with the associated mode shapes:

$$ \begin{align} \omega_0 &= \sqrt{\alpha} = 93.6278 & \Phi_0 &= [1, 1, 1, 1, 1, 1]^T \notag \\ \omega_1 &= \sqrt{\alpha+c} = 94.4168 & \Phi_1^c &= [1, 1, 0, -1, -1, 0]^T \notag \\ \omega_2 &= \sqrt{\alpha+3c} = 95.9753 & \Phi_1^s &= [1, \frac{1}{2}, -\frac{1}{2}, -1, -\frac{1}{2}, \frac{1}{2}]^T \notag \\ \omega_3 &= \sqrt{\alpha+4c} = 96.7451 & \Phi_2^c &= [1, -\frac{1}{2}, -\frac{1}{2}, 1, -\frac{1}{2}, -\frac{1}{2}]^T \notag \\ \omega_4 &= \sqrt{\alpha+5c} = 97.4511 & \Phi_3 &= [1, -1, 1, -1, 1, -1]^T \notag \end{align} \quad (41) $$

These mode shapes are depicted in Fig. 3.

## 4.2. Analytical study of nonlinear modes: free exact solutions

This section of the paper will focus on analytically determining the free solutions of system in Eq. (25). Results from the linear analysis will be used to determine solution forms for the nonlinear system in Eq. (36). Attention will be paid to the case of modes whose mode shapes contain just 0, 1 or −1 as a component (i.e. the $\Phi_p^c$ variables in Eq. (40)). For these mode shapes, all nonzero amplitude degrees of freedom actually vibrate with the same amplitude (in absolute value terms); consequently, it can be assumed that the adjustment of eigenfrequencies occurs naturally without having to change the modal shape. It will now be demonstrated that these linear modes are similar for the nonlinear system, meaning that their shapes do not vary with vibration amplitude. Let us suppose then that the solution vector $X$ assumes one of the following forms:

$$ \begin{align} \mathbf{X}_{(p=0)}(t) &= x_1(t)\mathbf{\Phi}_0 \quad \text{or} \quad \mathbf{X}_{(p=1)}(t) = x_1(t)\mathbf{\Phi}_1^c \notag \\ \mathbf{X}_{(p=2)}(t) &= x_1(t)\mathbf{\Phi}_2^c \quad \text{or} \quad \mathbf{X}_{(p=3)}(t) = x_1(t)\mathbf{\Phi}_3 \tag{42} \end{align} $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Mode shapes of the linear system: (a) mode $\Phi_0$; (b) mode $\Phi_1^c$; (c) mode $\Phi_1^s$; (d) mode $\Phi_2^c$; (e) mode $\Phi_2^s$; (f) mode $\Phi_3$.</figcaption></figure>
In this case, the system in Eq. (36) can be simplified into a single Duffing type equation:

$$x_1 + \omega_{0,p}^2 x_1 + \beta x_1^3 = 0 \quad (43)$$

with $\omega_{0,p}^2 = (\alpha + 2c(1 - \cos(\theta_p)))$ and $\theta_p = 2p\pi/6$ for $p = 0, 1, 2$ or $3$.

It can then be deduced that the linear mode shapes $\Phi_p^c$ for $p = 0, 1, 2$ or $3$ also serve as mode shapes for the nonlinear system. The nonlinear modes associated with these mode shapes are similar for the system in Eq. (36). Moreover, Eq. (43) enables identifying the relation between angular frequency and vibration amplitude for each of the four similar modes by means of Jacobi's elliptic functions (elliptic cosine and integrals) [29]. In the present case, the analytical expression of exact solutions to Eq. (43) is written as follows:

$$x_{1,p}(t) = A \operatorname{cn}(\lambda_p t; k_p) \quad (44)$$

where $cn(\cdot; \cdot)$ is the elliptic cosine, A the vibration amplitude and the quantities $\lambda_p, k_p$ are given by

$$\lambda_p^2 = \omega_{0,p}^2 + \beta A^2$$

$$k_p^2 = \frac{\beta A^2}{2(\omega_{0,p}^2 + \beta A^2)} \quad (45)$$

Lastly, the relation tying vibration amplitude A and vibration frequency $\omega_p$ is expressed by

$$\omega_p = \frac{\lambda_p \pi}{2K(k_p)} \quad (46)$$

where $K(k)$ is the complete elliptic integral of the first kind for the modulus $k$ as defined by

$$K(k) = \frac{2}{\pi} \int_{0}^{1} \frac{dx}{\sqrt{(1-x^2)(1-k^2x^2)}} \quad (47)$$

These analytical expressions will be used in order to validate our numerical results.

### 4.3. Numerical study of nonlinear modes for the 6 degree-of-freedom system

The previous section offered an effective display of free solutions from an analytical perspective. Attention will now be turned to a numerical calculation of the free solutions of Eq. (25).

#### 4.3.1. Convergence versus harmonic number

The effect of the selected harmonic numbers on solution convergence will first be studied by examining the zero-diameter mode. This step is equivalent to searching for solutions to Eq. (43) with a value of $p=1$. Three successive calculations were performed by retaining just 1, 2 or 3 odd harmonics. The results of this step are presented in Fig. 4a, where the vibration amplitude of the first degree of freedom has been plotted vs. frequency. It can be observed that the fact of considering the first and third harmonics is sufficient not only for solution convergence but also for an accurate approximation up until vibration amplitudes reaching an order of $h$. Moreover, the backbone curve obtained by HBM using two odd harmonics is compared with the exact solution obtained to Eq. (46). These results are presented in Fig. 4b which also reveals a strong level of agreement in results, as the relative error on frequency for a vibration amplitude equal to the thickness ($h=0.03$) is less than 0.1%.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Solution convergence and comparison of approximated solution with exact solution for the zero-diameter mode: (a) Convergence of the zero-diameter mode vs. harmonic number. (b) Comparison between the backbone approximated by HBM (two odd harmonics) and the exact solution.</figcaption></figure>
In light of these results, we will now be searching for the NNMs of Eq. (36) by assuming just the first and third harmonics, i.e. in the following form:

$$X(t) = A_1 \cos(\omega t) + A_3 \cos(3\omega t) \quad (48)$$

where $A_k$ are vectors of dimension 6.

In order to illustrate the backbone curves of this system, we will be using a representation within the energy-frequency plane analogous to that produced by Peeters in [3]. Let us note however a slight difference since the following relation will be introduced to express energy:

$$E(\mathbf{X}) = \sqrt{\sum_k \|\mathbf{A}_k\|} \quad (49)$$

which is directly expressed through the use of the HBM coefficients. This representation allows clearly displaying the minimum energy required for nonlinear effects to become significant, yet exhibits the disadvantage of concealing the individual behavior of degrees of freedom. To overcome this shortcoming, we will also provide the mode shape at some different frequencies.

### 4.3.2. Nonlinear normal modes

“Natural” NNMs: Nonlinear modes stemming from linear modes are depicted on the energy-frequency plane in Figs. 5 and 6. It is clearly observed that these modes are tangent to the linear normal modes and moreover that nonlinearity exerts a strong hardening effect (high frequency offset as energy increases) due to the retained Ritz function (Section 2).

For double modes (i.e. one and two diameters), it is observed that the backbone curves of both modes associated with the same frequency differ (Fig. 6) as a result of a number of sectors of the form $n = 2m$, with $m$ being an odd number.

It has been verified that modes associated with the forms $\Phi_i^c$ ($i = 0, 1, 2$, or 3) are similar; their mode shapes do not vary with vibration amplitude. On the other hand, the modes associated with forms $\Phi_i^s$ and $\Phi_i^b$ are nonsimilar, and the first

<figure><img src="image_6.png" /><figcaption>Fig. 5. Backbone curves and mode shapes for: (a) Zero-diameter NNM. (b) Three-diameter NNM.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Backbone curves and mode shapes for: (a) One-diameter NNM. (b) Two-diameter NNM.</figcaption></figure>
tends toward a mode shape where all amplitudes are equal, while the second tends toward a vibration mode that localizes motion on just two blades (Fig. 6). These results are in close agreement with the analytical study conducted in the previous section.

**Localized nonlinear mode:** In this section we concentrate on the NNM with three nodal diameters (NNM3D). By using the procedure described in Section 3.3, we detected several branching point bifurcations and we computed bifurcated branches. They are represented in Fig. 7 (dashed curves) along with the backbone curve of NNM3D (continuous curve). The mode shapes of bifurcated branches for a frequency about 190 rad s⁻¹ are depicted in Fig. 8. The bifurcated branches correspond to localized NNMs which means that only one (strongly localized) or several (weakly localized) sectors have a non-negligible amplitude of motion, the other sectors remaining virtually motionless [6]. It turns out that these bifurcated branches correspond to the (approximative) localized solution computed in Section 4.1. It can thus be concluded that solutions localized on 1, 2, 3 or 4 sectors originate as of the three-diameter mode through branching point bifurcations.

All these localized modes will prove useful when studying the system's forced vibrations. More specifically, it will be seen that during an excitation that assumes the three-diameter mode form, bifurcated solutions corresponding to a localized motion can be identified; moreover, these solutions are positioned around the backbones calculated in this section.

**Traveling waves:** Should the terms be kept in sine form as part of the development in Eq. (30), and by conveniently choosing the form of the solution at low vibration amplitude, it then becomes possible to calculate the NNMs that correspond to waves propagating in both the clockwise and counterclockwise directions. In the case of a six degree-of-freedom system, these nonlinear modes correspond to one and two-diameter rotating waves. In addition, in the case of a two-diameter traveling wave, a branching point bifurcation can be detected; this bifurcation gives rise to a traveling wave localized on three blades. The backbone curves of traveling waves in the energy-frequency plane are shown in Fig. 9; afterward, Fig. 10 depicts the temporal evolution of the various sectors, clearly indicating that the sectors vibrate with a constant phase shift between two consecutive sectors (the traveling waves shown here are propagating in a clockwise direction). Fig. 10c illustrates the special case of the traveling wave localized on three sectors.

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Backbone curve of the three-diameter mode (−) and its localized bifurcations (−−) (a: one sector, b: two sectors, c: three sectors, d: four sectors).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Mode shapes of three-diameter mode bifurcations for a frequency ω = 190 rad s⁻¹: (a) Localized mode on one sector. (b) Localized mode on two sectors. (c) Localized mode on three sectors. (d) Localized mode on four sectors.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Backbone curve of the nonlinear modes corresponding to traveling waves.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Temporal evolution of traveling waves: (a) 1-Diameter traveling wave. (b) 2-Diameter traveling wave. (c) 2-Diameter localized traveling wave.</figcaption></figure>
## 5. Forced response

Previous results reveal tremendous complexity in the free vibration behavior of such structures. It is obvious that this complexity will reappear in the case of the forced regime. We propose examining the impact of the various types of loadings (modal load, localized, detuned) on the structure of these responses; the focus thus lies on the forced response of the 6 degree-of-freedom system. The numerical parameters used to calculate the forced response are those defined in Eq. (28). With the aim of deriving curves of finite amplitudes, a damping term is added to Eq. (25) which then becomes

$$[M]\ddot{X} + [C]\dot{X} + [K]X + \beta X^3 = F(t) \quad (50)$$

with $[C] = \delta[I]$ and $\delta$ being a constant defined by $\delta = \omega_0/200$.

The right-hand side of the equation $F(t)$ is considered to be of the following form:

$$F(t) = A_F \cos(\Omega t) \quad (51)$$

where vector $A_F$ represents both the form and amplitude of the force.

The solutions to this problem are then calculated using both the harmonic balance method and continuation algorithm presented above. After convergence studies, we are assuming herein that determining the solutions with two odd harmonics is sufficient to provide us with accurate amplitude and stability results. The forced solutions of the system in Eq. (50) will thus be sought in the form

$$X(t) = A_1 \cos(\Omega t) + B_1 \sin(\Omega t) + A_2 \cos(3\Omega t) + B_2 \cos(3\Omega t) \quad (52)$$

Results will be depicted in frequency-energy plots, with the energy defined as follows:

$$E(\mathbf{X}) = \sqrt{\sum_k \|\mathbf{A}_k\|^2 + \sum_k \|\mathbf{B}_k\|^2} \quad (53)$$

### 5.1. Zero-diameter modal excitation

This section of the paper entails exciting the system along the zero-diameter mode, thus the excitation force can be written as

$$F(t) = a_f \Phi_0 \cos(\Omega t) \quad (54)$$

where $\Phi_0$ is the shape of the zero diameter mode (Eq. (41)).

The influence of force amplitude $a_f$ on both the form and stability of the response will be studied herein, results are presented in Fig. 11. It is noted that for small excitation amplitudes, the nonlinear response is analogous to the linear response and moreover that this solution is stable for all excitation frequencies. By increasing the force amplitude, a zone of instability is exposed between two turning points. By increasing the force amplitude even further, a second zone of instability is created after the second turning point (Fig. 12). This second zone originates via a branching point bifurcation. Through the use of the method described in Section 3, it is possible to calculate a secondary solution branch, which corresponds to a weakly localized vibration mode, where the sectors no longer vibrate in phase (Fig. 13). The bifurcated

<figure><img src="image_12.png" /><figcaption>Fig. 11. Nonlinear forced response for an excitation on the zero-diameter mode for different values of $a_f$ (—: stable, ···: unstable, - - -: linear response, --: backbone curve).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. Zoom of the bifurcated branches that arise when $a_f = 1/4$ (—: stable, ···: unstable, - - -: linear response, --: backbone curves).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Temporal evolution and amplitude at the bifurcated solution peak for an excitation on the zero-diameter mode with $a_f = 1/4$: (a) Temporal evolution. (b) Absolute value of amplitude.</figcaption></figure>

solution is stable over a few points, then loses its stability through a Hopf bifurcation, and finally a stable part is identified at the bifurcated solution "peak".

It is observed that even simple diameter excitations reveal dynamic complexity. Consequently, it may be assumed that localized excitations or detuned excitations are capable of exposing even greater complexity and varying types of behavior. The following section is intended to verify this assumption.

## 5.2. Localized excitation

This section adopts the stance that the excitation force is only acting on a single substructure (the first one) and thus assumes the following expression:

$$F(t) = a_f(1, 0, 0, 0, 0)^T \cos(\Omega t) \qquad (55)$$

The forced response is calculated for various values of amplitude $a_f$. These results are presented in Figs. 14 and 15. For purposes of comparison, the linear system response of Eq. (37) has also been depicted on these figures.

For $a_f \approx 2.5/16$, we can observe that the nonlinear system response is topologically equivalent to the linear system response. When $a_f \approx 1/4$, the appearance of a stable forced localization can be detected, meaning that for a conveniently chosen excitation frequency (between 98 and 100 rad s⁻¹), only the excited sector vibrates with a noticeable amplitude. Furthermore, this forced localization actually occurs around the backbone curve of the strongly localized mode obtained in the free vibration part, which is in agreement with what has already been demonstrated in [7].

When $a_f \approx 3.5/16$ (Fig. 15), it is possible to show a different form of the previous frequency response functions. These solutions are represented by closed curves that are completely dissociated from the nonlinear FRF and positioned around the backbone curves of the strongly localized mode. By studying the stability of this secondary solution, its potential
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Nonlinear forced response for a localized excitation on Sector 1 for different values of $a_f$ (—: stable, ···: unstable, - - -: linear response, --: backbone curves).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Zoom on the closed curve solution and the forced localization that arise in the case of a localized excitation (—: stable, ···: unstable, - - -: linear response, --: backbone curves).</figcaption></figure>

stability is apparent for certain excitation frequency values, and moreover the vibration amplitude (related to the energy) of this solution may be well above that of the basic nonlinear FRF.

This second solution, dissociated from the first, cannot be accessed by a single continuation; instead, to obtain the closed curve solution, an initiation point must be chosen on the localized FRF branch for a fixed force amplitude such as $a_f = \frac{1}{4}$, next, a sequential continuation with decreasing force amplitude is performed on this point with a fixed $\Omega$. The subsequent step consists of applying the arc-length continuation algorithm at the selected point to determine the solution in closed curve form (Fig. 15).

The identification of solutions in closed curve form leads to assuming that the forced localization phenomenon is due to the fusion of two distinct solution families (i.e. classical nonlinear response and closed curve), which coexist for intermediate excitation amplitude values and which overlap when $a_f$ is greater than a threshold value. Moreover, solutions in closed curve form may be stable for relatively high vibration amplitudes and for a seemingly weakly nonlinear FRF behavior. This remark is an important one since such solutions might be overlooked when applying a classical continuation approach, which could lead to erroneous structural design.

### 5.3. Detuned zero-diameter mode excitation

If, in the case of excitations, assuming a diameter mode form offers a "simpler" response than in the case of a localized excitation, it must be kept in mind that actual loadings are often imperfect and correspond more to perturbed diameter
excitations than to a perfect diameter. Our system will now be excited in a detuned manner, meaning that a perturbation is added to the excitation force. Let us consider that this perturbation only affects the force amplitude, yielding the following form:

$$F(t) = (\Phi + \epsilon)\cos(\Omega t) \qquad (56)$$

where $\epsilon$ is a vector characterizing the perturbation.

Due to the scalar products induced by the harmonic balance method, disturbing the force in this manner is equivalent to adding a perturbation to the $G(Y, \omega)$ function (Eq. (33)) similar to calculating the branching points. In this case however, perturbation is present throughout the duration of the calculation. The branching points will thus disappear and bifurcations will become a part of the solution, which can then be tracked by continuation. Furthermore, the fact of detuning the force makes it possible to evaluate the detuning level at which localization phenomena first appear. Our study focuses on the zero-diameter vibration mode, yet equivalent results can still be obtained for other vibration modes. We assume that just the first force component is in fact disturbed, a case that yields

$$F(t) = a_f(1 + \epsilon_1, 1, 1, 1, 1)^T \cos(\Omega t) \qquad (57)$$

where $\epsilon_1$ represents the force detuning percentage. Throughout the remaining discussion, it will be considered that $a_f = \frac{1}{4}$.

The simulations are launched with a 1% detuning level. It is easily observed in Fig. 16, that bifurcation of the $\epsilon_1 = 0$ case (Fig. 11) has since been included in the curve generated by continuation, as was expected.

By increasing the detuning percentage, amplitudes also rise and as of a 26% detuning level, it is possible to identify a solution uncorrelated with the primary solution, as in the case of localized excitation. The amplitudes have been

<figure><img src="image_17.png" /><figcaption>Fig. 16. Nonlinear forced response for an excitation on the detuned zero-diameter mode with several values of detuning $\epsilon_1$ (—: stable, ···: unstable, - - -: linear response, - -: backbone curves).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Zoom on the closed curve solution and the forced localization that arise in the case of a detuned zero diameter excitation (—: stable, ···: unstable, - - -: linear response, - -: backbone curves).</figcaption></figure>
represented in Figs. 16 and 17 for a 40% detuning level. This second family of solutions corresponds to a strongly localized vibration mode; moreover, these secondary solutions merge with the primary solution for a detuning level in the vicinity of 70%.

In this case, the closed curves appear for a detuning level of about 26%, which is quite significant. However, one has to keep in mind that the appearance of closed curves is closely related to the damping ratio: a smaller damping ratio would have cause the closed curves to appear at a smaller detuning level. Given that these closed curves contain stable zones corresponding to a substantial vibration amplitude (at least equal to the amplitude of the case $c_1 = 0$), they must be taken into account during structural design. A change in initial conditions or an operating shock could actually trigger such solutions and place a heavy load on one of the substructures.

## 6. Conclusion

In the research presented herein, a calculation was performed of both the free and forced responses of a structure with cyclic symmetry subjected to geometric nonlinearities. Setting up the system of equations, as discretized according to the Ritz method, led to a hardening nonlinear system. The nonlinear modes as well as the forced system responses were calculated using the harmonic balance method, while solution stability was determined by means of Floquet's theory. Besides the nonlinear modes originating from linear normal modes, other possible vibration modes were exposed; these stemmed from natural nonlinear modes via bifurcation. Some of these bifurcated modes correspond to somewhat localized vibration modes. Moreover, the presence of free solutions in the form of traveling waves has been underscored. These solutions are, in turn, capable of bifurcating into a form of localized traveling wave. In the case of a localized excitation, in addition to classical responses, a type of solution was identified in the form of closed curves. These closed curves lie adjacent to the classical solutions as the force amplitude increases, thereby inducing the phenomenon of forced localization. For an excitation on diameter modes, branching point type bifurcations were revealed, and these might also correspond to stable and localized vibration modes. Lastly, in the case of excitation detuned with respect to the zero-diameter mode, secondary solutions were discovered along with the forced localization phenomenon. This study has served to utilize the complexity of responses for structures with cyclic symmetry subjected to geometric nonlinearities. More specifically, we found closed curve solutions of potential importance during the design of these mechanical systems. An estimation of these solutions, obtained through a loading level analysis of the various nonlinear FRF, could be improved by conducting a systematic search for all solutions, e.g. by means of homotopy.

## References

[1] S. Samaranayake, Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry, *Journal of Sound and Vibration* 206 (1) (1997) 39-60.

[2] A.F. Vakakis, Dynamics of a nonlinear periodic structure with cyclic symmetry, *Acta Mechanica* 95 (1992) 197-226.

[3] F. Georgiades, M. Peeters, G. Kerschen, J. Golinval, Modal analysis of a nonlinear periodic structure with cyclic symmetry, *AIAA Journal* 47 (2009) 195-216.

[4] A.F. Vakakis, A multiple-scales analysis of nonlinear, localized modes in a cyclic periodic system, *Journal of Applied Mechanics* 60 (1993) 388-397.

[5] Y.J. Yana, P. Cuia, H. Hao, Vibration mechanism of a mistuned bladed-disk, *Journal of Sound and Vibration* 317 (2008) 294-307.

[6] A. Vakakis, *Normal Mode and Localization in Nonlinear Systems*, Wiley-Interscience, 1996.

[7] A.F. Vakakis, M. King, A very complicated structure of resonances in a nonlinear system with cyclic symmetry: non linear forced localization, *Nonlinear dynamics* 7 (1995) 85-104.

[8] O.O. Bendiksen, Localization phenomena in structural dynamics, *Chaos, Solitons and Fractals* 11 (2000) 1621-1660.

[9] M.E. King, A.F. Vakakis, Asymptotic analysis of nonlinear mode localization in a class of coupled continuous structures, *International Journal of Solids and Structures* 32 (1995) 1161-1177.

[10] R. Bhattacharyya, P. Jain, A. Nair, Normal mode localization for a two degree of freedom system with quadratic and cubic non linearities, *Journal of Sound and Vibration* 249(5) (2002) 909-919.

[11] M. Peeters, G. Kerschen, R. Viguié, G. Sérandour, J. Golinval, Nonlinear normal modes, Part ii: toward a practical computation using continuation technique, *Mechanical System and Signal Processing* 23 (2009) 195-216.

[12] R. Benamar, M. Bennouna, R. White, The effect of large vibration amplitudes on the mode shapes and natural frequencies of thin elastic structures, Part ii: fully clamped rectangular isotropic plates, *Journal of Sound and Vibration* 164 (2) (1993) 295-316.

[13] M. Amabili, Theory and experiments for large-amplitude vibrations of rectangular plates with geometric imperfections, *Journal of Sound and Vibration* 291 (2006) 539-565.

[14] M. Amabili, *Nonlinear Vibration and Stability of Shells and Plates*, Cambridge University Press, 2004.

[15] K. Liew, C. Wang, pb2-Rayleigh-Ritz method for general plate analysis, *Engineering structures* 15 (1) (1993) 55-60.

[16] D. Laxalde, F. Thouverez, J.J. Sinou, J.P. Lombard, Qualitative analysis of forced response of blisks with friction ring dampers, *European Journal of Mechanics and Solids* 36 (2007) 676-687.

[17] D. Laxalde, F. Thouverez, Complex nonlinear modal analysis for mechanical systems: application to turbomachinery bladings with friction interfaces, *Journal of Sound and Vibration* 322 (2009) 1009-1025.

[18] R. Lewandowski, Computational formulation for periodic vibration of geometrically nonlinear structures—Part 1: theoretical background, *International Journal of Solids Structures* 34 (15) (1997) 1925-1947.

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

[20] A. Nayfey, B. Balanchandran, *Applied Nonlinear Dynamics*, Wiley-Interscience, 1995.

[21] G. von Groll, D.J. Ewins, The harmonic balance method with arc length continuation in rotor stator contact problems, *Journal of Sound and Vibration* 241(2) (2000) 223-233.

[22] E. Sarrouy, Analyse globale de systèmes mécaniques non linéaires. application à la dynamique des rotors, PhD Thesis, Ecole Centrale Lyon, 2008.

[23] P. Ribeiro, M. Petyt, Non-linear free vibration of isotropic plates with internal resonance, *International Journal of Non-Linear Mechanics* 35 (2000) 263-278.
[24] A.F. Vakakis, Nonlinear normal mode and their application in vibration theory: an overview, *Mechanical System and Signal Processing* 11 (1) (1996) 3–22.

[25] G. Kerschen, M. Peeters, J. Golinval, A.F. Vakakis, Nonlinear normal modes, Part i: a useful framework for the structural dynamicist, *Mechanical System and Signal Processing* 23 (2009) 170–194.

[26] R. Rosenberg, On nonlinear vibration of systems with many degrees of freedom, *Advance in Applied Mechanics* (1966) 155–242.

[27] S. Shaw, C. Pierre, Normal modes for non-linear vibratory systems, *Journal of Sound and Vibration* 164 (1) (1993) 85–124.

[28] C. Touze, O. Thomas, A. Chaigne, Hardening/softening behaviour in non linear oscillation of structural systems using non linear normal modes, *Journal of Sound and Vibration* 273 (2004) 77–101.

[29] T. Okabe, T. Kondou, Improvement to the averaging method using the Jacobian elliptic function, *Journal of Sound and Vibration* 320 (2009) 339–364.