# Vibration Analysis of a Nonlinear System With Cyclic Symmetry

Aurélien Grolet, Fabrice Thouverez

► To cite this version:

Aurélien Grolet, Fabrice Thouverez. Vibration Analysis of a Nonlinear System With Cyclic Symmetry. Journal of Engineering for Gas Turbines and Power, 2011, 133 (2), pp.022502. 10.1115/1.4001989 . hal-01977265

HAL Id: hal-01977265

https://hal.science/hal-01977265

Submitted on 24 Jan 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, Fabrice Thouverez

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

This work is devoted to the study of nonlinear dynamics of structures with cyclic symmetry under geometrical nonlinearity using the harmonic balance method (HBM). In order to study the influence of the nonlinearity 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, 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. This 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

## 1 Introduction

The study considers both free and forced vibrations of structures with a cyclic symmetry under geometrical nonlinearity. This class of systems appears in the model of a bladed disk or a 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 pair 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, localization 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 led to many scientific papers [4,5]. Unlike the 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 nonlinearity 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 disappears. Traveling wave motions have also been detected in such systems [1,7]. In recent works, Peeters [1,8] used a shooting method coupled with a continuation algorithm to study the NNMs of a system with a cyclic symmetry. Similar and nonsimilar 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 disk using the harmonic balance method (HBM) coupled with an arc length continuation. This study emphasizes the numerous bifurcations that can happen in this kind of system. Attention has been paid to the localization phenomenon

and particularly to the link between nonsymmetric loading and localization. The effect of the force amplitude on the solutions is also studied.

## 2 Simplified Model

Consider a structure with a cyclic symmetry made of n identical substructures called sectors. As a result of its dimension and 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 et al. [9] and Amabili [10] focused on the case of thin rectangular plates with large deflection. The case of circular shells has been studied by Touzé et al. [11], and the case of beams has been studied by Lewandowski [12,13]. Here, the computation of our simplified mathematical model for a structure with cyclic symmetry is described. 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)).

Let us consider a single 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)$ ($(Oy)$, $(Oz)$) is denoted by $u$ ($v$, $w$) (Fig. 1(b)).

The Love–Kirchhoff hypothesis for the displacement, the use of the Von Karman nonlinear strain-displacement relationships, and the standard bidimensional Hooke law [14] lead to the following expression for the elastic strain energy U for one single plate:

$$ U = \frac{1}{2} \frac{Eh}{(1 - \nu^2)} \int_{0}^{L_x} \int_{-L_y/2}^{L_y/2} \left( \left( \frac{\partial w}{\partial x} \right)^2 \right)^2 + \left( \frac{\partial w}{\partial y} \right)^2 + 2\nu \left( \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2 \frac{1}{2} \left( \frac{\partial w}{\partial y} \right)^2 \right) + \frac{1 - \nu}{2} \left( \frac{\partial w}{\partial x} \frac{\partial w}{\partial y} \right)^2 dxdy + \frac{1}{2} \frac{Eh^3}{12(1 - \nu^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 + 2\nu \frac{\partial^2 w}{\partial x^2} \frac{\partial^2 w}{\partial y^2} $$
<figure><img src="image_2.png" /></figure>

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

$$ + \frac{1 - \nu}{2} \left( 2 \frac{\partial^2 w}{\partial x \partial y} \right)^2 dxdy \quad (1) $$

where $E$ is the Young modulus and $\nu$ is the Poisson ratio. By neglecting rotary inertia and the term $\dot{u}^2 + \dot{v}^2$ (since it is in $O(h^2)$, which is small compared with $\dot{w}^2$, which is in $O(1)$), 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{w}^2 dxdy \quad (2) $$

In the remainder of this paper, only a harmonic force, orthogonal to the plate, localized in $(x_f, y_f)=(L_x, 0)$ is considered. The work $W$ due to such an excitation is given by

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

where $\Omega$ is the excitation frequency, $t$ is the time, and $A_f$ is the force amplitude.

Let us 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$ is its strain energy, $T^j$ is its kinetic energy, and $W^j$ is 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 Eqs. (1)–(3). The coupling between substructures is modeled by massless linear stiffness of value $k$. This stiffness is positioned in a point $(x_r, y_r)=(L_x/4, 0)$ for all plates. The energy $V^j$ of such stiffness between 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 \text{for } 1 \le j \le n \quad (4) $$

with convention $j+1=1$ if $j=n$ (cyclic symmetry).

The total energies $U_t$, $T_t$, $V_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$ is 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 $w^j(x, y, t) = \sum_{i=1}^{N} \lambda_i^j(t) \Phi_i(x, y)$, 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}$.

In this paper, one proposes 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 = (x/L_x)^2$; it satisfies the boundary condition clamped at $x=0$. By using Lagrange's equations one finally obtained the following equations of motion:

$$ \ddot{X} + KX + \beta X^3 = F(t) \quad (5) $$

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

$$ 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 (6) $$

where $\alpha$, $c$, $\beta$, and $F$ are defined by

$$ \alpha = \frac{5}{3} \frac{Eh^2}{\rho L_x^4 (1 - \nu^2)}, \quad \beta = \frac{8E}{\rho L_x^4 (1 - \nu^2)}, \quad m = \frac{\rho h L_x L_y}{5}, \\ c = \frac{1}{256} \frac{k}{m}, \quad F^j(t) = \frac{A_j^j \cos(\omega t)}{m} \quad (7) $$

## 3 Resolution Method: Harmonic Balance Method

**3.1 Principle of the HBM.** The HBM is a widely spread method for solving nonlinear differential equations. One of the advantages of this method is that it can weakly or strongly treat a nonlinear system in the same way. It can be used for the problem with friction [16,17] or geometrical nonlinearity [12,18]. The HBM consists in finding periodic solutions $X$ of Eq. (5) of the form

$$ X(t) = A_0 + \sum_{k=1}^{N_h} A_k \cos(k\omega t) + B_k \sin(k\omega t) \quad (8) $$

where $N_h$ is the number of retained harmonics. Substituting Eq. (8) in Eq. (5) 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 $(A_k)_{0 \le k \le N_h}$, $(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. On one hand, the higher $N_h$ is, the better the solution. But on the other hand, when $N_h$ is too high, computations may require much time and memory. However, in many cases, few harmonics are needed to ensure good convergence leading to a reasonable system size.
**3.2 Numerical Procedure: Arc Length Continuation.** In a general way, the algebraic system obtained by the HBM can be reformulated as

$$G(Y, \omega) = 0 \qquad (9)$$

where $Y$ is a vector of dimension $n(2N+1)$ containing the unknowns $(A_k)_{0 \le k \le N}$ and $(B_k)_{1 \le k \le N}$, and $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 numerically finding a series of points $(Y^i, \omega^i)_{0 \le i \le M}$, which verifies the convergence criterion $G(Y^i, \omega^i) \le \epsilon$, where $\epsilon$ is a small parameter that determines the accuracy of the solution. In an arc length continuation, solutions are parametrized by the arc length $s$, so that $(Y^i, \omega^i) = (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 corrects the solution until the convergence criterion is achieved. In the study, a tangent predictor and the Newton-Raphson corrector have been chosen. Further explanations about continuation techniques can be found in Ref. [19].

**3.3 Stability of Solutions and Bifurcations.** The stability of periodic solutions is determined by the Floquet theory. It needs the computation of the monodromy matrix and of its eigenvalues. A detailed presentation of the method is given in Ref. [19]. Here, the focus is on the bifurcation that the system can withstand. There are two ways 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] = \partial G_l / \partial Y$ of the system in Eq. (9): 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 = \partial G_l / \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 are two ways 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 has been detected with the continuation algorithm, one removes the perturbation and carries 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 the system in Eq. (5). The natural approach is to use the NNMs, which have led to numerous publications. A survey about NNMs is given in Refs. [4,5]. The most used definition of NNM is the one given by Rosenberg in Ref. [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 such as 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 Eq. (5) will first be studied. This will give starting points for the study of the nonlinear system. Then, the search for the NNMs of the system in Eq. (5) by the HBM, coupled with an arc length continuation, will be carried out.

Numerical values for the physical parameters are introduced: $L_x=1.5$ m, $L_y=0.3$ m, $h=0.03$ m, $E=210$ GPa, $\nu=0.3$, and $k=8 \times 10^5$ N m$^{-1}$, which correspond to the following values for the parameters of Eq. (5): $\alpha=8.7662 \times 10^3$ s$^{-2}$, $c=148.36$ s$^{-2}$, and $\beta=4.6752 \times 10^7$ m$^{-2}$ s$^{-2}$.

**4.1.1 Modal Analysis of the Underlying Linear System.** The underlying linear system of problem in Eq. (5) is given by $\ddot{X} + KX = F(t)$. 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] \\ \omega_1 &= \sqrt{\alpha + c} = 94.42\Phi_1^c = [1, 1, 0, -1, -1, 0] \\ \Phi_1^s &= [1, 1/2, -1/2, -1, -1/2, 1/2] \\ \omega_2 &= \sqrt{\alpha + 3c} = 95.98\Phi_2^c = [1, -1, 0, 1, -1, 0] \\ \Phi_2^s &= [1, -1/2, -1/2, 1, -1/2, -1/2] \\ \omega_3 &= \sqrt{\alpha + 4c} = 96.75\Phi_3 = [1, -1, 1, -1, 1, -1] \end{aligned} \qquad (10) $$

**4.1.2 Nonlinear Normal Modes.** In this section, the nonlinear normal modes of the system in Eq. (5) are computed with the algorithm described earlier. NNM branches are computed by starting from the corresponding linear normal mode at low amplitudes. 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(X) = \sqrt{\|A_0\|^2 + \sum_k (\|A_k\|^2 + \|B_k\|^2)} \qquad (11)$$

where $A_k$ and $B_k$ are the HBM coefficients defined in Eq. (8). Since the nonlinearity is odd and because there is no damping in Eq. (5), only odd harmonics and cosine terms are retained in the development of $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 harmonics on the zero nodal diameter mode showed that a good approximation 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 the system in Eq. (5) are sought of the form $X(t) = A_1 \cos(\omega t) + A_3 \cos(3\omega t)$.

**4.1.3 Natural Nonlinear Normal Modes.** We call a natural NNM a NNM that is directly born from a LNM at low magnitudes of vibration. Since there are six LMNs, there will be six natural NNMs. They are represented in Figs. 2 and 3. 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 non-similar NNMs. NNMs born from the LNM shapes $\Phi_i^c$ ($i=0, 1, 2$, or 3) defined in Eq. (10) are similar NNMs, and NNMs born from the LNM shapes $\Phi_i^s$ ($i=1, 2$) are nonsimilar NNMs. For double modes, the two backbone curves associated with the $p$ nodal diameter mode ($p=1$ or 2) are different (Fig. 3). This is because the number of sectors is of the form $n=2m$ with $m$ odd. As we shall see, other NNM branches bifurcate from these backbone curves.
<figure><img src="image_4.png" /><figcaption>Fig. 2 Backbone curves and deformed shapes for natural NNMs with (a) 0 or (b) three nodal diameters</figcaption></figure>

### 4.1.4 Bifurcated NNMs.
This section focuses on the NNM with three nodal diameters (NNM3D). By using the procedure described in Sec. 3.3, several branching point bifurcations have been detected and bifurcated branches have been computed. They are represented in Fig. 4 (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 also depicted in Fig. 4. 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 localiza-

tion phenomenon was already observed in a linear mistuned sys- tem [22], but here, and in general nonlinear system with cyclic symmetry, it occurs without a structural disorder. The spatial con- finement of the energy causes the responses of some sectors to be high and might lead to a premature failure of the blades.

### 4.1.5 Traveling Waves.
All the NNMs computed so far corre- spond to standing wave motions in the sense that the DOFs vi- brate in a synchronous way as in Rosenberg’s definition of NNMs. Some papers [1,7] refer to traveling wave motions in which vibra- tions are not synchronous anymore. For Vakakis [7], these travel-

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

<figure><img src="image_6.png" /><figcaption>Fig. 4 Backbone curve of NNM3D (–) and its localized bifurcations (––): (a) one sector, (b) two sectors, (c) three sectors, and (d) four sectors</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5 Backbone curves of traveling wave motions and time series of the two nodal diameters localized traveling wave</figcaption></figure>

ing 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 NNM computation. Sine terms are re-introduced in the HBM development of solution X in order to take into account the phase difference. Therefore, X is now sought of the following form:

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

Well chosen starting points enable the computation of traveling waves propagating in the clockwise or anticlockwise direction. The backbone curves of traveling waves are given in Fig. 5(a). By monitoring the Jacobian matrix, as indicated in Sec. 3.3, a branch that bifurcates from the traveling wave with two nodal diameters has been found. 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 a non-negligible amplitude. To illustrate the phase difference between coordinates, we plotted in Fig. 5(b) the time series of the six blades for the two nodal diameter, localized traveling wave motions.

## 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,2,7]. 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 out for the structure with six DOFs, which correspond to Eq. (5). Numerical parameters

are the same as in the previous section. In order to have a finite amplitude response, a damping term is added to Eq. (5), which becomes

$$\ddot{X} + \delta\dot{X} + KX + \beta X^3 = F(t) \quad (13)$$

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

The right hand side of Eq. (13) is considered to be of the form $F(t)=A_F \cos(\Omega t)$, where $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. The stability of solutions is investigated using the Floquet theory, and bifurcations are detected as described in Sec. 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 Eq. (13) are sought of the form $X(t)=A_1 \cos(\Omega t)+B_1 \sin(\Omega t) +A_3 \cos(3\Omega t)+B_3 \cos(3\Omega t)$. 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 a 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 the system in Eq. (13) 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

$$F(t) = A_F \cos(\Omega t) \quad \text{with} \quad A_F = a_f \Phi_0 \quad (14)$$

The effect of the force amplitude $a_f$ on the form and on the stability of solutions is studied. For low force amplitudes (Fig.

<figure><img src="image_8.png" /><figcaption>Fig. 6 Nonlinear forced response for a low engine order with zero nodal diameter (–, stable; ..., unstable; –.–, linear response)</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7 Nonlinear forced response for a low engine order with zero nodal diameter with $a_f=1/4$ (–, stable; ..., unstable; –.–, linear response)</figcaption></figure>

6(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. 6(b)), an unstable zone is generated between two turning points. When the force amplitude is increased again (Fig. 7), 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 Sec. 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 corresponds to a weakly localized motion (framed areas of Fig. 7).

## 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[1,0,0,0,0,0]^T \cos(\Omega t) \quad (15)$$

The forced response has been computed for several values of the amplitude $a_f$. Results are depicted in Figs. 8 and 9(a). Because of symmetry, only the response for four of the six sectors has been represented. For values of amplitudes from $a_f=0$ to $a_f=a_f^s$

$\approx 0.8/4$, the nonlinear response is topologically equivalent to the response of the underlying linear system (Fig. 8(a)). When $a_f$ is larger than $a_f^{loc} \approx 0.9125/4$, stable forced localization occurs (Fig. 9(a)): For a well chosen excitation frequency (between 98 rad s⁻¹ and 100 rad s⁻¹), only the excited sector vibrates 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 has been computed in the previous section (free vibration). These results are in agreement with those of Vakakis presented in Ref. [23]. When $a_f^s \le a_f \le a_f^{loc}$ (Fig. 8(b)), it is possible to find another type of solution. These new solutions are represented by a 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 much larger amplitudes 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 a fixed frequency and with the force amplitude $a_f$ as a parameter has been used. The starting point was taken from the

<figure><img src="image_10.png" /><figcaption>Fig. 8 Nonlinear forced response for a localized excitation on sector 1 (–, stable; ..., unstable; –.–, linear response)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9 Nonlinear forced response for a localized excitation on sector 1 (–, stable; ..., unstable; –.–, linear response) and principle of the closed curve detection</figcaption></figure>

part of the curve that corresponds 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. The method for closed curve detection is illustrated in Fig. 9(b). The new kind of solution makes us think that forced localization is the result of the fusion of the closed curves 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.

### 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 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 (16)$$

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 Eq. (9). 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[1 + \varepsilon_1, 1, 1, 1, 1]^T \cos(\Omega t) \quad (17)$$

where $\varepsilon_1$ corresponds to the detuning percentage of the force with a zero nodal diameter. In the remainder of this paper, the force amplitude $a_f$ will be set to $a_f=1/4$. Simulations start with $\varepsilon=1\%$. Results are represented in Fig. 10. As expected, the bifur-

<figure><img src="image_12.png" /><figcaption>Fig. 10 Nonlinear forced response for a detuned excitation with $\varepsilon_1=1\%$ (–, stable; ..., unstable; –.–, linear response)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11 Nonlinear forced response for a detuned excitation (–, stable; ..., unstable; –, linear response)</figcaption></figure>

cated branch of the case ε=0% (Fig. 7) is now a part of the solution. When ε₁ becomes larger than 26%, it is possible to compute a secondary solution that corresponds to a closed curve. This kind of solution is represented in Fig. 11(a) for ε₁=40%. This closed curve has stable parts, and it merges with the basic nonlinear solution for ε₁=70% (Fig. 11(b)). In this example, the closed curve has been detected when the detuning is larger than 26%, which is a quite significant detuning level. However, the threshold for the detection of the 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 has been obtained by a Rayleigh-Ritz method. Solutions were computed by the harmonic balance method, stability was investigated using the 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 a closed curve in the amplitude frequency diagram, and it merges with the basic 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 dynamics of a system with cyclic symmetry under geometrical 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.

## Nomenclature

$A_k, B_k$ = HBM coefficients

$F(t)$ = time domain external forces

$G(Y, \omega)$ = algebraic system given by the HBM

$J_y, J_\omega$ = Jacobian of G with respect to $\dot{Y}$ or $\omega$

$K_{nl}$ = nonlinear stiffness

$M, K$ = mass and linear stiffness matrices

$X(t)$ = time domain vector of DOF

$Y$ = vector of HBM coefficient

$L_x, L_y, h$ = plate length, width, and thickness

$T, U$ = kinetic and strain energy

$V$ = strain energy of a linear stiffness

$W$ = work due to external forces

$u, v, w$ = displacements

$\lambda_i^j$ = Ritz coordinates (DOF)

$\Phi_i^c, \Phi_i^s$ = linear deformed shapes

$\Phi_i$ = Ritz shape functions

$\rho$ = density

$\omega, \Omega$ = free frequency, excitation frequency

## References

[1] Georgiades, F., Peeters, M., Kerschen, G., and Golinval, J., 2009, “Modal Analysis of a Nonlinear Periodic Structure With Cyclic Symmetry,” AIAA J., 47, pp. 1014–1025.

[2] Vakakis, A. F., 1993, “A Multiple-Scales Analysis of Nonlinear, Localized Modes in a Cyclic Periodic System,” ASME J. Appl. Mech., 60, pp. 388–397.

[3] Samaranayake, S., 1997, “Subharmonic Oscillations in Harmonically Excited Mechanical Systems With Cyclic Symmetry,” J. Sound Vib., 206(1), pp. 39–60.

[4] Vakakis, A. F., 1997, “Non-Linear Normal Modes (NNMs) and Their Applications in Vibration Theory: An Overview,” Mech. Syst. Signal Process., 11(1), pp. 3–22.

[5] Kerschen, G., Peeters, M., Golinval, J., and Vakakis, A. F., 2009, “Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist,” Mech. Syst. Signal Process., 23, pp. 170–194.

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

[7] Vakakis, A. F., 1992, “Dynamics of a Nonlinear Periodic Structure With Cyclic Symmetry,” Acta Mech., 95, pp. 197–226.

[8] Peeters, M., Kerschen, G., Vigui, R., Srandour, G., and Golinval, J., 2009, “Nonlinear Normal Modes, Part II: Toward a Practical Computation Using Continuation Technique,” Mech. Syst. Signal Process., 23, pp. 195–216.

[9] Benamar, R., Bennouna, M., and White, R., 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,” J. Sound Vib., 164(2), pp. 295–316.

[10] Amabili, M., 2006, “Theory and Experiments for Large-Amplitude Vibrations of Rectangular Plates With Geometric Imperfections,” J. Sound Vib., 291, pp. 539–565.

[11] Touzé, C., Thomas, O., and Chaigne, A., 2004, “Hardening/Softening Behaviour in Non Linear Oscillation of Structural Systems Using Non Linear Normal Modes,” J. Sound Vib., 273, pp. 77–101.

[12] Lewandowski, R., 1997, “Computational Formulation for Periodic Vibration of Geometrically Nonlinear Structures—Part I: Theoretical Background,” Int. J. Solids Struct., 34(15), pp. 1925–1947.

[13] Lewandowski, R., 1997, “Computational Formulation for Periodic Vibration of
Geometrically Nonlinear Structures—Part 2: Numerical Strategy and Examples," Int. J. Solids Struct., 34(15), pp. 1949–1964.
[14] Ribeiro, P., and Petyt, M., 2000, "Non-Linear Free Vibration of Isotropic Plates With Internal Resonance," Int. J. Non-Linear Mech., 35, pp. 263–278.
[15] Liew, K., and Wang, C., 1993, "pb-2 Rayleigh-Ritz Method for General Plate Analysis," Eng. Struct., 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," Eur. J. Mech. A/Solids, 26, pp. 676–687.
[17] Laxalde, D., and Thouverez, F., 2009, "Complex Nonlinear Modal Analysis for Mechanical Systems: Application to Turbomachinery Bladings With Friction Interfaces," J. Sound Vib., 322, pp. 1009–1025.
[18] Ribeiro, P., and Petyt, M., 1999, "Nonlinear Vibration of Plates by the Hier-

archical Finite Element and Continuation Methods," Int. J. Mech. Sci., 41, pp. 437–459.
[19] Nayfeh, A. H., and Balachandran, B., 1995, *Applied Nonlinear Dynamics*, Wiley-Interscience, New York.
[20] Rosenberg, R., 1966, "On Nonlinear Vibration of Systems With Many Degrees of Freedom," Adv. Appl. Mech., 9, pp. 155–242.
[21] Shaw, S., and Pierre, C., 1993, "Normal Modes for Non-Linear Vibratory Systems," J. Sound Vib., 164(1), pp. 85–124.
[22] Yan, Y. J., Cui, P. L., and Hao, H. N., 2008, "Vibration Mechanism of a Mistuned Bladed-Disk," J. Sound Vib., 317, pp. 294–307.
[23] Vakakis, A. F., and King, M., 1995, "A Very Complicated Structure of Resonances in a Nonlinear System With Cyclic Symmetry: Non Linear Forced Localization," Nonlinear Dyn., 7, pp. 85–104.