# Vibration Analysis of a Nonlinear System With Cyclic Symmetry

Aurelien Grolet, Fabrice Thouverez

► To cite this version:

Aurelien Grolet, Fabrice Thouverez. Vibration Analysis of a Nonlinear System With Cyclic Symmetry. ASME Turbo Expo 2010: Power for Land, Sea, and Air, Jun 2010, Glasgow, United Kingdom. pp.917-929, 10.1115/GT2010-22681. hal-02424290

HAL Id: hal-02424290

https://hal.science/hal-02424290

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
# VIBRATION ANALYSIS OF A NONLINEAR SYSTEM WITH CYCLIC SYMMETRY

Aurélien Grolet<sup>a,*</sup>, Fabrice Thouverez<sup>a</sup>

(a) École Centrale de Lyon, Laboratoire de Tribologie et Dynamique des Systèmes
36, avenue Guy de Collongue, 69134 Ecully Cedex, France

E-mail: Aurelien.Grolet@ec12009.ec-lyon.fr

## ABSTRACT

This work is devoted to the study of non linear dynamics of structures with cyclic symmetry under geometrical nonlinearity using the harmonic balance method (HBM). In order to study the influence of the non-linearity due to large deflection of blades a simplified model has been developed. It leads to nonlinear differential equations of the second order, linearly coupled, in which the nonlinearity appears by cubic terms. Periodic solutions in both free and forced cases are sought by the HBM coupled with an arc length continuation and stability analysis.

In this study, a specific attention has been paid to the evaluation of nonlinear modes and to the influence of excitation on dynamic responses. Indeed, several cases of excitation have been analyzed: punctual one and tuned or detuned low engine order. The paper shows that for a localized, or sufficiently detuned, excitation, several solutions can coexist, some of them being represented by closed curves in the Frequency-Amplitude domain. Those different kinds of solution meet up when increasing the force amplitude, leading to forced nonlinear localization. As the closed curves are not tied with the basic nonlinear solution they are easily missed. They were calculated using a sequential continuation with the force amplitude as a parameter.

**Keywords:** Geometrical nonlinearity, nonlinear normal modes, bifurcation, localization.

## NOMENCLATURE

A<sub>k</sub>, B<sub>k</sub> HBM coefficients.

F(t) Time domain external forces.

G(Y, ω) Algebraic system given by the HBM.

K<sub>nl</sub> Nonlinear stiffness.

M, K Mass and Linear stiffness matrices.

T, U Kinetic and strain energy.

V Strain energy of a linear stiffness.

W Work due to external forces.

X(t) Time domain vector of dof.

Y Vector of HBM coefficient.

L<sub>x</sub>, L<sub>y</sub>, h Plate length, width and thickness

J<sub>y</sub>, J<sub>ω</sub> Jacobian of G with respect to Y or ω.

u, v, w Displacements.

Φ<sub>i</sub><sup>c</sup>, Φ<sub>i</sub><sup>s</sup> Linear deformed shapes.

Φ<sub>i</sub> Ritz shape functions.

λ<sub>i</sub><sup>j</sup> Ritz coordinates (dof).

ρ Density.

ω, Ω Free frequency, Excitation frequency.

## 1 Introduction

The study considers both free and forced vibrations of structures with cyclic symmetry under geometrical nonlinearity. This class of systems appears in the model of bladed disc or space antenna [1], [2]. It leads to nonlinear differential equations of the second order, linearly coupled, in which the nonlinearity appears by cubic terms.

In the linear case, most of the natural frequencies appear in pairs due to the perfect symmetry of the problem. These natural frequencies are related to deformed shapes with nodal diameters [3]. For weakly coupled and weakly mistuned systems, lo-
calization can take place leading to motions that are confined on only a few substructures.

In the nonlinear case, the study of free vibrations relies on the definitions of nonlinear normal modes (NNMs) which have leads to many scientific papers [4] [5]. Unlike linear system, the number of NNMs can exceed the number of dofs. The other NNMs arise from bifurcations. In both free and forced cases, if the ratio between the coupling of the substructure and the non-linearity is small, Vakakis [6] showed that nonlinear localization can take place in a perfectly symmetric system. When this ratio increases, bifurcations occur and nonlinear localization disappear. Traveling waves motions have also been detected in such systems [7], [1].

In recent works, Peeters [8] [1] used a shooting method coupled with a continuation algorithm to study the NNMs of a system with cyclic symmetry. Similar and non similar NNMs have been found. Moreover, he studied the modal interactions between modes. He showed that these interactions can occur even if the natural frequencies of the modes are not commensurable, and he detected a countable infinity of such interactions.

The aim of this paper is to study both free and forced nonlinear vibrations of a bladed disc, using the HBM coupled with an arc length continuation. This study emphasizes the numerous bifurcation that can happen in this kind of system. Attention has been paid to the localization phenomenon, and particularly to the link between non symmetric loading and localization. The effect of the force amplitude on the solutions is also studied.

## 2 Simplified model

Consider a structure with cyclic symmetry made of n identical substructures called sectors. As a result of its dimension, its materials, and because of the external effort, such sectors can experiment large deflections leading to geometrical nonlinearity.

Many papers have studied the effect of geometrical nonlinearity on thin structures. Among others, Benamar [9] and Amabili [10], focused on the case of thin rectangular plates with large deflection. The case of circular shells have been studied by Touzé [11], and the case of beams by Lewandowski [12] [13].

Here the computation of our simplified mathematical model for a structure with cyclic symmetry is describe. In this paper, each sector is modeled by a thin rectangular plate clamped on one edge to a fixed frame. The coupling between the substructure is realized by a linear stiffness (fig.1 a). First, only one sector will be considered in order to obtain the expression of the different kinds of energy. These preliminary results will then be applied to the system with cyclic symmetry.

### 2.1 Preliminary considerations for one sector

This section focuses on a single sector. The retained hypothesis will be presented and expressions for the different kinds of

<figure><img src="image_2.png" /><figcaption>Figure 1. Retained model and system of coordinates for a rectangular plate</figcaption></figure>

energy will be given.

Consider a plate P with dimensions $L_x, L_y$, and thickness h in a Cartesian system of coordinates $(O, x, y, z)$, clamped on the edge $(x = 0)$. The displacement of a point with coordinates $(x, y, z)$ in the direction $(Ox)$ (respectively $(Oy), (Oz)$) is denoted by $u$ (respectively $v, w$) (fig.1 b).

**Elastic strain energy of a plate** The Love-Kirchhoff hypothesis for the displacement are used so that the displacement $u$ is expressed as:

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

In addition, the use of the Von Karman nonlinear strain-displacement relationships and the standard bidimensional Hooke's law [14], leads to the following expression for the elastic strain energy U:
$$ U = \frac{1}{2} \frac{Eh}{(1 - v^2)} \int_{0}^{L_x} \int_{-L_y/2}^{L_y/2} \left( \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2 \right)^2 + \left( \frac{1}{2} \left( \frac{\partial w}{\partial y} \right)^2 \right)^2 + 2v \left( \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2 \right) \frac{1}{2} \left( \frac{\partial w}{\partial y} \right)^2 + \frac{1 - v}{2} \left( \frac{\partial w}{\partial x} \frac{\partial w}{\partial y} \right)^2 dxdy \\ + \frac{1}{2} \frac{Eh^3}{12(1 - v^2)} \int_{0}^{L_x} \int_{-L_y/2}^{L_y/2} \left( \frac{\partial^2 w}{\partial x^2} \right)^2 + \left( \frac{\partial^2 w}{\partial y^2} \right)^2 + 2v \frac{\partial^2 w}{\partial x^2} \frac{\partial^2 w}{\partial y^2} + \frac{1 - v}{2} \left( 2 \frac{\partial^2 w}{\partial x \partial y} \right)^2 dxdy \quad (2) $$

where E is the Young's modulus and v is the Poisson's ratio.

**Kinetic energy and external load** By neglecting rotary inertia, the kinetic energy T of a rectangular plate is given by:

$$ T = \frac{1}{2} \rho h \int_{x=0}^{L_x} \int_{y=-L_y/2}^{L_y/2} (\dot{u}^2 + \dot{v}^2 + \dot{w}^2) dxdy \quad (3) $$

In this paper the terms $\dot{u}^2 + \dot{v}^2$ will be omitted since it is in $O(h^2)$ which is small compared to $\dot{w}^2$ which is in $O(1)$. The kinetic energy is therefore given by:

$$ T = \frac{1}{2} \rho h \int_{x=0}^{L_x} \int_{y=-L_y/2}^{L_y/2} \dot{w}^2 dxdy \quad (4) $$

The work W due to external forces is written as:

$$ W = \int_{x=0}^{L_x} \int_{y=-L_y/2}^{L_y/2} (f_x u + f_y v + f_z w) dxdy \quad (5) $$

where $f_x, f_y, f_z$ represent the distributed forces per unit area acting respectively in (Ox), (Oy), (Oz) directions. In this study only a localized harmonic force orthogonal to the plate is considered, thus $f_x = f_y = 0$. The external load $f_z$ of a localized force is given by:

$$ f_z = A_f \cos(\Omega t) \delta(x - x_f) \delta(y - y_f) \quad (6) $$

where $\Omega$ is the excitation frequency, $t$ is the time, $\delta$ is the Dirac function, $A_f$ is the force amplitude and $(x_f, y_f) = (L_x, 0)$. In this case, equation (5) can be rewritten as:

$$ W = w(L_x, 0) A_f \cos(\Omega t) \quad (7) $$

## 2.2 Equation of motion for the structure with cyclic symmetry

Let's come back to the structure with cyclic symmetry. For $1 \le j \le n$, $w^j$ denotes the transverse displacement of the plate number j, $U^j$ its strain energy, $T^j$ its kinetic energy, and $W^j$ the work due to external forces on this plate. $U^j, T^j$ and $W^j$ are obtained by substituting $w$ by $w^j$ in equations (2), (4), (7). The coupling between substructures is modeled by massless linear stiffness of value $k$. These stiffness are positioned in a point $(x_r, y_r) = (L_x/4, 0)$ for all plates. The energy $V^j$ of such stiffness between the plates j and j + 1 is given by:

$$ V^j = \frac{1}{2} k (w^j(x_r, y_r) - w^{j+1}(x_r, y_r))^2 \quad (8) $$
for $1 \le j \le n$ with convention $j+1 = 1$ if $j = n$

The total energies $U_t, U_t, U_t$ and $W_t$ are then given by the sum over the number of plates $n$ of the different local energies $U^j, T^j, V^j$ and $W^j$.

The discretization of the transverse displacement $w^j$ in now introduced. This discretization is done by a Rayleigh-Ritz method [15], and it allows us to have a simplified model of the system. Transverse displacements $w^j$ are interpolated by the following expression:

$$ w^j(x, y, t) = \sum_{i=1}^{N} \lambda_i^j(t) \Phi_i(x, y) \quad (9) $$

where $(\Phi_i)$ are kinematically admissible shape functions, $\lambda_i^j$ is the contribution of the shape function $\Phi_i$ in the displacement $w^j$, and $N$ is the number of shape functions used for the interpolation. This method leads to a discretized problem with $nN$ unknowns $(\lambda_i^j)_{1 \le i \le N, 1 \le j \le n}$. Equations of motion are then given by Lagrange's equations (Appendix).

In this paper, we propose to study a system with six identical substructures ($n=6$), of which displacements are interpolated by a single Ritz shape function ($N=1$). This will lead to a problem with six degrees of freedom which will be used for numerical computation.

The retained shape function for the example is $\Phi = (\frac{x}{L_x})^2$, it satisfies the boundary condition clamped at $x=0$. By using Lagrange's equations we finally obtained the following equations of motion:

$$ \dot{\mathbf{X}} + \mathbf{K}\mathbf{X} + \beta\mathbf{X}^3 = \mathbf{F}(t) \quad (10) $$

with $\mathbf{X} = (\lambda_j)_{1 \le j \le 6}$, the convention $\mathbf{X}^3 = (x_j^3)_{1 \le j \le 6}$, and $\mathbf{K}$ a
matrix defined by:

$$ \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} \quad (11) $$

where $\alpha$, $c$, $\beta$, and $\mathbf{F}$ are given in appendix

# 3 Resolution method: Harmonic Balance Method

## 3.1 Principle of the HBM

The HBM is a widely spread method for solving non linear differential equations. One of the advantages of this method is that it can treat weakly or strongly non linear system in the same way. It can be used for problem with friction [16] [17], or geometrical nonlinearity [12], [18].

Consider the following system of non linear differential equations

$$ \mathbf{M}\ddot{\mathbf{X}} + \mathbf{K}\mathbf{X} + \mathbf{K}_{nl}(\mathbf{X}) = \mathbf{F}(t) \quad (12) $$

where $\mathbf{X}$ is a vector of unknowns of dimension $n$, $\mathbf{M}$ and $\mathbf{K}$ the mass and linear stiffness matrix of the system, $\mathbf{K}_{nl}$ the non linear forces, and $\mathbf{F}$ the external forces.

The HBM consists in finding of periodic solutions $\mathbf{X}$ of the form:

$$ \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) \quad (13) $$

where $N_h$ is the number of retained harmonics.

Substituting equation (13) in equation (12), and projecting the result on the trigonometric base $(1, (\cos(k\omega t), \sin(k\omega t)))_{1\le k\le N_h}$, one obtains a set of $n(2N_h + 1)$ algebraic nonlinear equations with $n(2N_h + 1) + 1$ unknowns $(\mathbf{A}_k)_{0\le k\le N_h}, (\mathbf{B}_k)_{1\le k\le N_h}$ and $\omega$.

The most important parameter in this method is the number $N_h$ of retained harmonics. This number is not known a priori, so convergence studies are needed to ensure a good estimation of the solution. In one hand the higher $N_h$ is the better the solution. But in another hand, when $N_h$ is too high, computations can require many time and memory. However, in many cases, few harmonics are needed to ensure a good convergence leading to reasonable system size.

## 3.2 Numerical procedure: Arc length continuation

In a general way, the algebraic system obtains by the HBM can be reformulated as following:

$$ \mathbf{G}(\mathbf{Y}, \omega) = 0 \quad (14) $$

where $\mathbf{Y}$ is a vector of dimension $n(2N + 1)$ containing the unknown $(\mathbf{A}_k)_{0\le k\le N}$ and $(\mathbf{B}_k)_{1\le k\le N}$, and $\mathbf{G}$ is a nonlinear function taking its values in a space of dimension $n(2N + 1)$. This system is going to be solved by an arc length continuation method. Continuation methods consist in finding numerically a series of points $(\mathbf{Y}^i, \omega^i)_{0\le i\le M}$ which verifies the following criterion :

$$ \mathbf{G}(\mathbf{Y}^i, \omega^i) \le \varepsilon \quad (15) $$

where $\varepsilon$ is a small parameter which determines the accuracy of the solution. In an arc-length continuation, solutions are parameterized by the arc-length $s$, so that $(\mathbf{Y}^i, \omega^i) = (\mathbf{Y}(s_i), \omega(s_i))$. The method takes place in two steps : a predictor step which gives an estimation of the solution and a corrector step which correct the solution until the criterion (15) is achieved. In the study a tangent predictor and the Newton-Raphson corrector have been chosen. Further explanations can be found in [19]. A sum up of the continuation algorithm is given in Fig.2. This algorithm has been implemented in the Matlab software.

<figure><img src="image_3.png" /><figcaption>Figure 2. Algorithm for the computation of solutions by Arc-Length continuation</figcaption></figure>

## 3.3 Stability of solutions and bifurcations

The stability of periodic solutions is determined by Floquet theory. It needs the computation of the monodromy matrix and
of its eigenvalues. A detailed presentation of the method is given in [19].

Here, the focus is on the bifurcation that the system can withstand. There is two way of detecting a bifurcation. First, by monitoring the eigenvalues of the monodromy matrix: if some eigenvalues go outside the unit circle then there is a bifurcation, the type of the bifurcation is determined by the way that the eigenvalues leave the unit circle. The second method to detect bifurcations is by monitoring the determinant $|J_y|$ of the jacobian matrix $[J_y] = \frac{\partial G}{\partial Y}$ of system (14): when this determinant is zero at some point, then the matrix $[J_y]$ is singular and this point is a bifurcation point. The type of bifurcation is determined by the range of the matrix $[J_y J_\omega]$ (where $J_\omega = \frac{\partial G}{\partial \omega}$). For a turning point, the range of $[J_y J_\omega]$ is $n$, and there is no particular treatment to apply since the arc length continuation can handle turning points. For a branching point, the range of $[J_y J_\omega]$ is at most $n-1$, so there is at least two tangent vectors $T = [T_y T_\omega]^t$ such that $[J_y J_\omega]T = 0$. There is two way of computing the bifurcated branches. The first method is to add a small perturbation to the function $G(Y, \omega)$ in the vicinity of the bifurcation point. Since branching points are structurally unstable, the perturbation breaks the bifurcation. Once another branch have been detected with the continuation algorithm, one remove the perturbation and carry on the continuation [19]. The other method is to compute the eigenvectors $\Psi_i$ of the matrix $[J_y]$ associated with the zero eigenvalue at the bifurcation point. The eigenvectors $\Psi_i$ indicate the directions to be followed [14].

## 4 Free vibrations of a system with six dofs

The first objective of this paper is to estimate the free solutions of system (10). The natural approach is to use the nonlinear normal modes (NNMs) which have lead to numerous publications. A survey about NNMs is given in [4] and [5]. The most used definition of NNM is the one given by Rosenberg in [20], in which he defines the NNMs by the vibrations in unison of a free and undamped system. Other approaches to define NNMs exist, they are based on geometrical arguments like the invariant manifold definition proposed by Shaw and Pierre [21], or based on the normal form theory [11].

In this section, the underlying linear system of problem (10) will first be studied. This will give starting points for the study of the nonlinear system. Then, the search for the NNMs of system (10) by the HBM coupled with an arc length continuation will be carried out.

Numerical values for the physical parameters are introduced:

$$ 
\begin{aligned}
L_x &= 1,5m, \quad L_y = 0,3m, \quad h = 0,03m, \\
E &= 210\text{GPa}, \quad \nu = 0.3, \quad k = 8.10^5\text{N.m}^{-1}
\end{aligned} 
\quad (16) 
$$

which correspond to the following values for the parameters of

equation (10):

$$ \alpha = 8,7662 \cdot 10^3 s^{-2}, \quad c = 148,36 s^{-2}, \quad \beta = 4,6752 \cdot 10^7 m^{-2} \cdot s^{-2} \quad (17) $$

**Modal analysis of the underlying linear system**
The underlying linear system of problem (10) is given by:

$$ \ddot{X} + KX = F(t) \quad (18) $$

Linear normal modes (LNMs) are given by the resolution of the eigenvalue problem $K\Phi = \omega^2\Phi$. Matrix $K$ have single and double eigenvalues, but double eigenvalues represent the majority [3]. With each double eigenvalue are associated two distinct deformed shapes which are not uniquely defined. The mode shapes are characterized by their number of nodal diameters $p$. The eigenvalues and eigenvectors of matrix $K$ are given by (the linear frequencies are given in rad.s$^{-1}$):

$$ 
\begin{aligned}
\omega_0 &= \sqrt{\alpha} = 93.63 & \Phi_0 &= [1, 1, 1, 1, 1, 1] \\
\omega_1 &= \sqrt{\alpha+c} = 94.42 & \Phi_1^c &= [1, 1, 0, -1, -1, 0] \\
& & \Phi_1^s &= [1, \frac{1}{2}, -\frac{1}{2}, -1, -\frac{1}{2}, \frac{1}{2}] \\
\omega_2 &= \sqrt{\alpha+3c} = 95.98 & \Phi_2^c &= [1, -1, 0, 1, -1, 0] \\
& & \Phi_2^s &= [1, -\frac{1}{2}, -\frac{1}{2}, 1, -\frac{1}{2}, -\frac{1}{2}] \\
\omega_3 &= \sqrt{\alpha+4c} = 96.75 & \Phi_3 &= [1, -1, 1, -1, 1, -1]
\end{aligned} 
\quad (19) 
$$

The deformed shapes of equation (19) are depicted in Fig.3

**Nonlinear normal modes** In this section, the nonlinear normal modes of system (10) are computed with the algorithm described earlier. NNMs branches are computed by starting from the corresponding linear normal mode at low amplitude. These branches are termed the backbone curves, and they are represented in an Energy-Frequency plot as in [1]. The energy considered here is defined by:

$$ E(\mathbf{X}) = \sqrt{\|A_0\|^2 + \sum_k (\|A_k\|^2 + \|B_k\|^2)} \quad (20) $$

where $A_k$ and $B_k$ are the HBM coefficients defined in relation (13). Since the non linearity is odd and because there is no damping in system (10), only odd harmonics and cosine terms are retained in the development of $\mathbf{X}$ (the last condition corresponds to a phase condition in which all initial velocities are set to zero). A convergence study with the number of harmonic on the zero nodal diameter mode showed that a good approximation
<figure><img src="image_4.png" /><figcaption>Figure 3. Mode shapes of the linear system (18)</figcaption></figure>

can be obtained when retaining only the first and the third harmonic. This approximation holds for amplitudes of vibration up to the plate thickness $h$. We suppose that this remark holds for all NNMs. Therefore the NNMs of system (10) are sought of the form:

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

**Natural nonlinear normal modes** We call a natural NNM, a NNM which is directly born from a LNM at low magnitude of vibration. Since there are six LMNs, there will be six natural NNMs. They are represented in Fig.4 and Fig.5. The first noticeable feature in these figures is the Frequency-Energy (or amplitude) dependence of the NNMs. In addition to the dependence of their oscillation frequency, the NNMs may also have their modal shapes varying with the amplitude of vibration: these are nonsimilar NNMs. NNMs born from the LNM shapes $\Phi_i^c$ ($i = 0, 1, 2$, or 3) defined in equation (19) are similar NNMs, and NNMs born from the LNM shapes $\Phi_i^s$ ($i = 1, 2$) are nonsimilar NNMs.

For double modes, the two backbones curves associated with the $p$ nodal diameter mode ($p = 1$ or 2) are different (fig.5). This is because the number of sector is of the form $n = 2m$ with $m$ odd.

As we shall see, other NNM branches bifurcate from these backbone curves.

**Bifurcated NNMs** This section focus on the NNM with 3 nodal diameters (NNM3D). By using the procedure described in section 3.3, several branching point bifurcations have been detected and bifurcated branches have been computed. They are

<figure><img src="image_5.png" /><figcaption>Figure 4. Backbone curves and deformed shapes for natural NNMs with 0 (a) or 3 (b) nodal diameters</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Backbone curves and deformed shapes for natural NNMs with 1 (a) and 2 (b) nodal diameters</figcaption></figure>

represented in Fig.6 (dashed curves) along with the backbone curve of NNM3D (continuous curve). The mode shapes of bifurcated branches for a frequency about $190\text{rad.s}^{-1}$ are depicted in Fig.7. The bifurcated branches correspond to localized NNMs. It means that only one (strongly localized) or several (weakly localized) sectors have a non negligible amplitude of motion, the others remaining virtually motionless [6]. Such localization phenomenon was already observed in a linear mistuned system [22], but here, and in general nonlinear system with cyclic symmetry, it occurs without structural disorder. The spatial confinement of the energy causes the responses of some sectors to be high and might lead to premature failure of the blades.

**Traveling waves** All the NNMs computed so far correspond to standing wave motions in the sense that the dofs vibrate in a synchronous way as in Rosenberg's definition of NNMs. Some papers [7], [1] refer to traveling wave motions in which vibrations are not synchronous anymore. For Vakakis [7], these traveling waves arise from a 1:1 internal resonance between the two modes associated with the same number of nodal diameter.

Because of the phase difference between the dofs, a different phase condition is considered for traveling wave NNMs computation. Sine terms are re-introduce in the HBM development of solution $X$ in order to take into account the phase difference.
<figure><img src="image_7.png" /><figcaption>Figure 6. Backbone curve of NNM3D (−) and its localized bifurcations (−·−)(a: 1 sector, b: 2 sectors, c: 3 sectors, d: 4 sectors)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7. Deformed shapes of the bifurcated modes for a frequency ω = 190rad.s⁻¹</figcaption></figure>

Therefore **X** is now sought of the following form:

$$ \mathbf{X}(t) = \mathbf{A}_1 \cos(\omega t) + \mathbf{B}_1 \sin(\omega t) + \mathbf{A}_2 \cos(3\omega t) + \mathbf{B}_2 \sin(3\omega t) \quad (22) $$

The starting point, at low amplitude of vibration, for the continuation algorithm is given by:

$$ \begin{aligned} \mathbf{A}_1 &= \varepsilon \cos(\mathbf{\Phi}_p) & \mathbf{A}_2 &= \mathbf{0} \\ \mathbf{B}_1 &= \varepsilon \sin(-\mathbf{\Phi}_p) & \mathbf{B}_2 &= \mathbf{0} \end{aligned} \quad (23) $$

where *p* correspond to the number of nodal diameter (*p* = 1, or 2), ε is the amplitude of vibration (considered small) and

$\mathbf{\Phi}_p$ is a vector which coordinates (($\mathbf{\Phi}_p$)$_i$)$_{1\le i\le 6}$ are defined by:

$$ (\mathbf{\Phi}_p)_i = (i-1) \frac{2p\pi}{6} \quad (24) $$

These starting points enable the computation of traveling waves propagating in the anticlockwise direction. Their companions, propagating in the clockwise direction, can be obtained by replacing $\mathbf{\Phi}_p$ by $-\mathbf{\Phi}_p$ in equation 23.

The backbone curves of traveling waves are given in Fig.8 (curves a and b). By monitoring the jacobian matrix as indicated in section 3.3, a branch that bifurcate from the traveling wave with 2 nodal diameters has been found (fig.8 curve c). This bifurcated branch is quite noticeable because it can be seen as a weakly localized traveling wave since only three of the six sectors vibrate with non negligible amplitude. To illustrate the phase difference between coordinates, we plotted in fig.9 the time series of the six blades for the 2 nodal diameters localized traveling wave motions.

<figure><img src="image_9.png" /><figcaption>Figure 8. Backbone curves of traveling wave motions with one nodal diameter (a), two nodal diameters (b), two nodal diameters localized on three sectors (c)</figcaption></figure>

## 5 Forced vibrations

Results from the previous section (free vibrations) are in agreement with the literature regarding structure with cyclic symmetry under geometric (cubic) nonlinearity [1], [7] [2]. They highlight the complexity of the free dynamic of such structures. Obviously, this complexity is going to come back in forced vibrations. We here propose to study the impact of the different external loading on the forced response. Simulations are carried
<figure><img src="image_10.png" /><figcaption>Figure 9. Time series of the two nodal diameters localized traveling wave</figcaption></figure>

out for the structure with six dofs which correspond to system (10). Numerical parameters are the same as in the previous section. In order to have finite amplitude response, a damping term is added to system (10) which becomes:

$$ \ddot{\mathbf{X}} + \delta\dot{\mathbf{X}} + \mathbf{K}\mathbf{X} + \beta\mathbf{X}^3 = \mathbf{F}(t) \quad (25) $$

with $\delta$ a damping coefficient defined by $\delta = \omega_0/200$.

The right hand side of equation (25) is considered to be of the following form:

$$ \mathbf{F}(t) = \mathbf{A}_F \cos(\Omega t) \quad (26) $$

where $\mathbf{A}_F$ stands for the shape and magnitude of the external force and $\Omega$ stands for the excitation frequency.

Solutions are computed by the HBM. Stability of solutions is investigated using Floquet theory, and bifurcations are detected as described in section 3.3. Since the nonlinearity is odd, only odd harmonics are retained in the HBM approximation. A convergence study of the stability with the number of retained harmonics showed that keeping the first and the third harmonic is sufficient enough to ensure valid results on stability. Therefore, solutions of equation (25) are sought of the form:

$$ \mathbf{X}(t) = \mathbf{A}_1 \cos(\Omega t) + \mathbf{B}_1 \sin(\Omega t) + \mathbf{A}_2 \cos(3\Omega t) + \mathbf{B}_2 \cos(3\Omega t) \quad (27) $$

Three kinds of excitation will be considered. First, low engine order excitation, then localized excitation, and finally, detuned excitation. In all cases, a parametric study with the force

amplitude as parameter has been carried out. Because of symmetry, only one (or four) of the six sectors are going to be represented.

## 5.1 Low engine order excitation with zero nodal diameters

In this section, the forced response of system (25) is studied for an excitation $F(t)$ which takes the shape of the zero nodal diameters LNM. This means that $F(t)$ is of the form:

$$ \mathbf{F}(t) = \mathbf{A}_F \cos(\Omega t) \quad \text{with} \quad \mathbf{A}_F = a_f \mathbf{\Phi} \quad (28) $$

where $\mathbf{\Phi}$ is equal to $\mathbf{\Phi}_0$ of equation (19)

The effect of the force amplitude $a_f$ on the form and on the stability of solutions is studied. For low force amplitude (Fig.10 a), the response of the nonlinear system is close to the one of the underlying linear system, and the solution is stable for all frequencies $\Omega$. By increasing the force amplitude (Fig.10 b), an unstable zone is generated between two turning points. When the force amplitude is increased again (Fig.11), a second zone of instability is generated after the second turning point. This second zone takes birth through a branching point bifurcation. By using the branch switching method described in section 3.3, a bifurcated branch of solution has been computed. The bifurcated branch is stable for a few points at the beginning of the curve and then an instability is generated through a Hopf bifurcation. The area near the peak of the bifurcated solution is stable and correspond to a weakly localized motion (framed areas of Fig.11)

<figure><img src="image_11.png" /><figcaption>Figure 10. Nonlinear forced response for a low engine order with zero nodal diameter (—: stable, ···: unstable, - -: linear response)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 11. Nonlinear forced response for a low engine order with zero nodal diameter with $a_f = \frac{1}{4}$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>

## 5.2 Localized excitation

In this section, the external loading acts on only one sector (the first one), therefore the force is given by:

$$F(t) = a_f \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{pmatrix} \cos(\Omega t) \quad (29)$$

The forced response has been computed for several values of the amplitude $a_f$. Results are depicted from Fig.12 to Fig.15. Because of symmetry only the response for four of the six sectors has been represented.

For values of amplitude from $a_f = 0$ to $a_f = a_f^s \approx \frac{0.8}{4}$, the nonlinear response is topologically equivalent to the response of the underlying linear system (fig.12).

When $a_f$ is larger than $a_f^{loc} \approx \frac{0.9125}{4}$, stable forced localization occurs (fig.15): for a well chosen excitation frequency (between 98 and 100 rad.s⁻¹), only the excited sector vibrate with a non negligible amplitude, and this kind of motion is stable. Moreover the part of the curve corresponding to forced localization is positioned around the backbone curve of the strongly localized NNM which have been computed in the previous section (free vibration). These results are in agreement with those of Vakakis presented in [23].

When $a_f^s \le a_f \le a_f^{loc}$ (Fig.14 and Fig.13), it is possible to

find another type of solution. These new solutions are represented by closed curve in the Amplitude-Frequency diagram. They are not tied up with the basic nonlinear response, but they are positioned around the backbone curve of the strongly localized NNM. Stability analysis shows that some parts of the closed curves can be stable for points that have very larger amplitude than the basic nonlinear solution.

Since the closed curves are not tied up with the basic nonlinear response, they are unreachable with a unique continuation scheme. To compute the closed curves, a sequential continuation with fixed frequency and with the force amplitude $a_f$ as a parameter has been used. The starting point was taken from the part of the curve that correspond to forced localization (for $a_f \ge a_f^{loc}$), then the force amplitude was decreased by continuation until $a_f \le a_f^{loc}$. The results obtained this way have been used as starting points for the arc length continuation. This method is illustrated in Fig.16.

The new kind of solution makes us think that forced localization is the result of the fusion of the closed curve with the basic nonlinear solution when $a_f = a_f^{loc}$. Whereas the basic nonlinear solution seems to be weakly nonlinear since it is close to the linear solution, there exist stable solutions (the closed curve) which have a large vibration amplitude and which could be attained depending on initial conditions. This remark is important because the closed curves could be easily missed by a classical continuation scheme, thus leading to a wrong design of the structure.

<figure><img src="image_13.png" /><figcaption>Figure 12. Nonlinear forced response for a localized excitation on sector 1 with $a_f = \frac{2.5}{16}$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 13. Nonlinear forced response for a localized excitation on sector 1 with $a_f = \frac{3,222}{16}$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Nonlinear forced response for a localized excitation on sector 1 with $a_f = \frac{3,4}{16}$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>

## 5.3 Detuned low engine order with zero nodal diameter excitation

It has been seen that in the case of low engine order excitation the system response is simpler than in the case of localized excitation. However we have to keep in mind that realistic excitations are not perfectly symmetric. So now, the external forces are considered detuned, that is a perturbation is added to the low

<figure><img src="image_16.png" /><figcaption>Figure 15. Nonlinear forced response for a localized excitation on sector 1 with $a_f = \frac{1}{4}$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Illustration of the method used to compute closed curve solutions</figcaption></figure>

engine order excitation. The perturbation is assumed to act only on the force amplitude. The external forces are then given by:

$$F(t) = (\Phi + \varepsilon) \cos(\Omega t) \quad (30)$$

where $\Phi$ is a linear deformed shape and $\varepsilon$ is a vector characterizing the perturbation

Because of the projection in the HBM, adding a perturbation to the force is equivalent to adding a perturbation to the function $G(Y, \omega)$ defined in equation (14). Then the branching point
bifurcations which are structurally unstable will disappear [19]. Thus, bifurcated branches will be part of the nonlinear solution and will be obtained by the continuation method.

Here, only a low engine order with zero nodal diameters is considered, but similar results could be obtained with the other vibration modes. Finally, the perturbation is assumed to act only on the first coordinate of the force. Therefore we have:

$$F(t) = a_f \begin{pmatrix} 1 + \varepsilon_1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} \cos(\Omega t) \quad (31)$$

where $\varepsilon_1$ correspond to the detuning percentage of the force with zero nodal diameter. In the remainder of the paper the force amplitude $a_f$ will be set to $a_f = \frac{1}{4}$.

Simulations start with $\varepsilon = 1\%$. Results are represented in Fig.17. As expected the bifurcated branch of the case $\varepsilon = 0\%$ (Fig.11) is now a part of the solution.

When $\varepsilon_1$ becomes larger than 26% it is possible to compute a secondary solution which correspond to a closed curve. This kind of solution is represented in Fig.18 for $\varepsilon_1 = 40\%$. This closed curve have stable parts, and it merges to the basic nonlinear solution for $\varepsilon_1 = 70\%$ (Fig.19).

In this example the closed curve have been detected when the detuning is larger than 26% which is a quite significant detuning level. However, the threshold for the detection of closed curve is highly dependent on the damping of the system: for a smaller damping the closed curves would have appeared sooner.

## 6 Conclusion

In this study, the free and forced response of a system with cyclic symmetry under geometric nonlinearity has been computed. The discretized equation of has been obtained by a Rayleigh-Ritz method. Solutions were computed by the harmonic balance method, stability was investigated using Floquet theory and bifurcations were computed through a branch switching algorithm.

In the free case, in addition to natural nonlinear normal modes, other vibration modes were detected by studying the branching point bifurcation of the system. Some of the bifurcated solutions correspond to localized nonlinear modes. Traveling wave motions have also been detected in this study.

In the forced case, when the excitation is localized to one of the substructure or sufficiently detuned, in addition to the basic nonlinear response, another kind of solution has been detected. This secondary solution is represented by closed curve in the Amplitude Frequency diagram, and it merge with the basic

<figure><img src="image_18.png" /><figcaption>Figure 17. Nonlinear forced response for a detuned excitation with $\varepsilon_1 = 1\%$ (—: stable, ···: unstable, - -: linear response)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 18. Nonlinear forced response for a detuned excitation with $\varepsilon_1 = 40\%$ (—: stable, ···: unstable, - -: linear response)</figcaption></figure>

nonlinear solution when increasing the force amplitude, leading to forced nonlinear localization. For low engine excitations, we computed several bifurcated branches which also correspond to stable localized motions.

This study is based on a simple model, therefore the results remain at a theoretical level. However, this paper highlights the complex dynamic of system with cyclic symmetry under geomet-
<figure><img src="image_20.png" /><figcaption>Figure 19. Nonlinear forced response for a detuned excitation with $\epsilon_1 = 70\%$ (—: stable, ···: unstable, - - -: linear response)</figcaption></figure>

rical nonlinearity. Of importance is the detection of closed curve solutions which can be stable, and which may play an important role during the design of such a mechanical system.

## Appendix: Expression of parameters for the system with six dofs

Here, we give the definition of parameters of equation (10).

With the definition of the Ritz shape function $\Phi$ given in section 2.2, the transverse displacement is given by:

$$w^j(x,y,t) = \lambda^j(\frac{x}{L_x})^2 \quad \text{for } 1 \le j \le 6 \quad (32)$$

Lagrange's equations are given by:

$$\frac{d}{dt} \frac{\partial T_t}{\partial \lambda_i^j} + \frac{\partial U_t}{\partial \lambda_i^j} + \frac{\partial V_t}{\partial \lambda_i^j} = \frac{\partial W_t}{\partial \lambda_i^j} \quad \text{for } 1 \le i \le N, 1 \le j \le n \quad (33)$$

Then the different terms of Lagrange's equations (33) become:

$$\begin{aligned} \frac{d}{dt} \frac{\partial T}{\partial \lambda_j^i} &= m \ddot{\lambda}_j^i \\ \frac{\partial U}{\partial \lambda_j^i} &= m \alpha \lambda_j^i + m \beta (\lambda_j^i)^3 \\ \frac{\partial V}{\partial \lambda_j^i} &= m c (2 \lambda_j^i - \lambda^{j+1} - \lambda^{j-1}) \\ \frac{\partial W}{\partial \lambda_j^i} &= m F^j(t) \end{aligned} \quad (34)$$

With:

$$\begin{aligned} m &= \rho \int_P \Phi^2 dP = \frac{\rho h L_x L_y}{5} \\ \alpha &= \frac{5}{3} \frac{E h^2}{\rho L_x^4 (1-v^2)} \\ \beta &= \frac{\delta E}{\rho L_x^4 (1-v^2)} \\ c &= \frac{k}{m} \Phi(x_r, y_r)^2 = \frac{1}{256} \frac{k}{m} \\ F^j(t) &= \Phi(x_f^j, y_f^j) \frac{f^j(t)}{m} = \frac{f^j(t)}{m} \end{aligned} \quad (35)$$

## References

[1] F.Georgiades, M. Peeters, G. K. J. G., 2009. "Modal analysis of a nonlinear periodic structure with cyclic symmetry". *AIAA journal*, **47**, pp. 195–216.

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

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

[4] Vakakis, A. F., 1996. "Nonlinear normal mode and their application in vibration theory: An overview". *Mechanical system and signal processing*, **11**(1), pp. 3–22.

[5] G. Kerschen, M. Peeters, J. G. A. F. V., 2009. "Nonlinear normal modes, part i: A useful framework for the structural dynamicist". *Mechanical system and signal processing*, **23**, pp. 170–194.

[6] Vakakis, A., 1996. *Normal mode and localization in nonlinear systems*. Wiley-Interscience.

[7] Vakakis, A. F., 1992. "Dynamics of a nonlinear periodic structure with cyclic symmetry". *Acta Mechanica*, **95**, pp. 197–226.

[8] M. Peeters, G. Kerschen, R. V. G. J. G., 2009. "Nonlinear normal modes, part ii: toward a practical computation using continuation technique". *Mechanical system and signal processing*, **23**, pp. 195–216.

[9] R. Benamar, M.M.K. Bennouna, R. W., 1993. "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), pp. 295–316.

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

[11] C. Touze, O.Thomas, A., 2004. "Hardening/softening behaviour in non linear oscillation of structural systems using non linear normal modes". *Journal of sound and vibration*, **273**, pp. 77–101.

[12] Lewandowski, R., 1997. "Computational formulation for periodic vibration of geometrically nonlinear structures-
part 1:theoretical background”. *International Journal of Solids Structures*, **34**(15), pp. 1925–1947.

[13] Lewandowski, R., 1997. "Computational formulation for periodic vibration of geometrically nonlinear structures-part 2: numerical strategy and examples" *International Journal of Solids Structures*, **34**(15), pp. 1949–1964.

[14] P. Ribeiro, M. P., 2000. "Non-linear free vibration of isotropic plates with internal resonance". *International Journal of Non-Linear Mechanics*, **35**, pp. 263–278.

[15] K.M. Liew, C. W., 1993. "pb2-rayleigh-ritz method for general plate analysis". *Engineering structures*, **15**(1), pp. 55–60.

[16] Laxalde, D., Thouverez, F., Sinou, J. J., and Lombard, J. P., 2007. "Qualitative analysis of forced response of blisks with friction ring dampers". *European journal of mechanics and solids*, **36**, pp. 676–687.

[17] Laxalde, D., and Thouverez, F., 2009. "Complex nonlinear modal analysis for mechanical systems: application to turbomachinery bladings with friction interfaces". *Journal of sound and vibration*, **322**, pp. 1009–1025.

[18] P. Ribiero, M. P., 1999. "Nonlinear vibration of plates by the hierarchical finite element and continuation method". *International journal of mechanical science*, **41**, pp. 437–459.

[19] A.H. Nayfey, B. B., 1995. *Applied nonlinear dynamics*. Wiley-Interscience.

[20] Rosenberg, R., 1966. "On nonlinear vibration of systems with many degrees of freedom". *Advance in applied mechanics*, pp. 155–242.

[21] S.W. Shaw, C. P., 1993. "Normal modes for non-linear vibratory systems". *Journal of sound and vibration*, **164**(1), pp. 85–124.

[22] Y.J.Yana, P.L. Cuia, H. H., 2008. "Vibration mechanism of a mistuned bladed-disk". *Journal of Sound and Vibration*, **317**, pp. 294–307.

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