Vibration Analysis of a Nonlinear System With Cyclic Symmetry

This work is devoted to the study of nonlinear dynamics of structures with cyclic symme-try under geometrical nonlinearity using the harmonic balance method (HBM). In order to study the inﬂuence of the nonlinearity due to large deﬂection of blades, a simpliﬁed 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 solu-tions in both free and forced cases are sought by the HBM coupled with an arc length continuation and stability analysis. In this study, speciﬁc attention has been paid to the evaluation of nonlinear modes and to the inﬂuence 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 sufﬁciently 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.


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.

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: This work is devoted to the study of nonlinear dynamics of structures with cyclic symme-try 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 solu-tions 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 where E is the Young modulus and is the Poisson ratio. By neglecting rotary inertia and the term u 2 + v 2 ͑since it is in O͑h 2 ͒, which is small compared with ẇ 2 , which is in O͑1͒͒, the kinetic energy T of a rectangular plate is given by 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 where ⍀ 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 Յ j Յ 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 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 where ͑⌽ i ͒ are kinematically admissible shape functions, i j is the contribution of the shape function ⌽ 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 ͑ i j ͒ 1ՅiՅN,1ՅjՅ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 ⌽ = ͑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: with X = ͑ j ͒ 1ՅjՅ6 , the convention X 3 = ͑ j 3 ͒ 1ՅjՅ6 , and K a matrix defined by where ␣, c, ␤, and F are defined by 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 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͑kt͒ , sin͑kt͒͒ 1ՅkՅN h ͒, one obtains a set of n͑2N h +1͒ algebraic nonlinear equations with n͑2N h +1͒ + 1 unknowns ͑A k ͒ 0ՅkՅN h , ͑B k ͒ 1ՅkՅN h , and . 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.
(a) Simplified model of a bladed disc (b) Rectangular plate and system of coordinates 3.2 Numerical Procedure: Arc Length Continuation. In a general way, the algebraic system obtained by the HBM can be reformulated as where Y is a vector of dimension n͑2N +1͒ containing the unknowns ͑A k ͒ 0ՅkՅN and ͑B k ͒ 1ՅkՅ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 where 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 , i ͒ = ͑Y͑s i ͒ , ͑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 achieve. 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͔.

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 ͔ = ‫ץ‬GÕ‫ץ‬ 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 ͔͑ where J = ‫ץ‬GÕ‫ץ‬ ͒. For a turning point, the range of ͓J y J ͔ 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 ͔ is at most n − 1, so there is at least two tangent vectors T = ͓T y T ͔ t such that ͓J y J ͔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 , ͒ 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 ⌿ i of the matrix ͓J y ͔ associated with the zero eigenvalue at the bifurcation point. The eigenvectors ⌿ i indicate the directions to be followed ͓14͔.

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.

Modal Analysis of the Underlying Linear
System. The underlying linear system of problem in Eq. ͑5͒ is given by Ẍ + KX = F͑t͒. Linear normal modes ͑LNMs͒ are given by the resolution of the eigenvalue problem K⌽ = 2 ⌽. 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 ͒ 0 = ͱ ␣ = 93.63⌽ 0 = ͓1,1,1,1,1,1͔

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 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.W e 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͑t͒ + A 3 cos͑3t͒.

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 nonsimilar NNMs. NNMs born from the LNM shapes ⌽ i c ͑i =0,1 ,2 , or 3͒ defined in Eq. ͑10͒ are similar NNMs, and NNMs born from the LNM shapes ⌽ 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.

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 −1 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 system ͓22͔, but here, and in general nonlinear system with cyclic symmetry, it occurs without a structural disorder. The spatial confinement of the energy causes the responses of some sectors to be high and might lead to a 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 ͓1,7͔ refer to traveling wave motions in which vibrations are not synchronous anymore. For Vakakis ͓7͔, these travel-  Fig. 4 Backbone curve of NNM3D "-… and its localized bifurcations "-.-…: "a… one sector, "b… two sectors, "c… three sectors, and "d… four sectors 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: 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.

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 param-eters 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 The right hand side of Eq. ͑13͒ is considered to be of the form F͑t͒ = A F cos͑⍀t͒, where A F stands for the shape and magnitude of the external force and ⍀ 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͑⍀t͒ + B 1 sin͑⍀t͒ + A 3 cos͑3⍀t͒ + B 3 cos͑3⍀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.

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 The effect of the force amplitude a f on the form and on the stability of solutions is studied. For low force amplitudes ͑Fig. 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 ⍀. By increasing the force amplitude ͑Fig. 6͑b͒͒,a n 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͒.

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͑⍀t͒͑ 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 Ϸ 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 Ϸ 0.9125/ 4, stable forced localization occurs ͑Fig. 9͑a͒͒: For a well chosen excitation frequency ͑between 98 rad s −1 and 100 rad s −1 ͒, only the excited sector vibrates with a nonnegligible 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 Յ a f Յ 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 part of the curve that corresponds to forced localization ͑for a f Ն a f loc ͒; then, the force amplitude was decreased by continuation until a f Յ 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.

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 where ⌽ is a linear deformed shape and 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 , ͒ 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+ 1 ,1,1,1,1,1͔ T cos͑⍀t͒͑ 17͒ where 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 = 1%. Results are represented in Fig. 10. As expected, the bifur-  cated branch of the case =0% ͑Fig. 7͒ is now a part of the solution. When 1 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 1 = 40%. This closed curve has stable parts, and it merges with the basic nonlinear solution for 1 = 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.

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 , ͒ ϭ algebraic system given by the HBM J y , J ϭ Jacobian of G with respect to Y or 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 i j ϭ Ritz coordinates ͑DOF͒ ⌽ i c , ⌽ i s ϭ linear deformed shapes ⌽ i ϭ Ritz shape functions ϭ density , ⍀ϭfree frequency, excitation frequency