E Sarrouy

F Thouverez
email: fabrice.thouverez@ec-lyon.fr

Global search of non-linear systems periodic solutions: a rotordynamics application

Keywords: Global analysis, Non-linear, Homotopy, Polynomial, Rotordynamics, Squeeze-film damper

Introducing non-linearities into models contributes towards a better reality description but leads to systems having multiple solutions. It is then legitimate to look for all the solutions of such systems, that is to have a global analysis approach. However no effective method can be found in literature for systems described by more than two or three degrees of freedom. We propose in this paper a way to find all T -periodic solutions -where T is known -of a nonlinear dynamical system. This method is compared to three other approaches and is shown to be the most efficient on a Duffing oscillator. As a more complex example, a rotor model including a squeeze-film damper is studied and a second branch of solutions is exhibited.

Introduction

Global analysis of dynamical systems is of great interest in industrial development as well as in research area.

Introduction of non-linear organs in dynamical modeling is one of the ways to a better reality description and so to a better behavior prediction; though it induces multiple solutions of varied nature (periodic, quasi-periodic...). One need to get them all in order to dimension a structure properly. Then, when trying to get the solutions of a set of nonlinear dynamical equations, two approaches, local or global, are possible. Local methods mostly consist in choosing a startpoint and having many correction steps until a solution point is found. On the contrary, global analysis methods consider the whole state space (or a bounded part of it) and try to find all the solutions contained in it.

Many methods aim at global analysis of non-linear systems of equations. They all share the same problem: their computational efficiency. In this paper a new approach is proposed for locating all T -periodic solutions. It is based on the transformation of the dynamical system into a non-differential one followed by a polynomial approximation of the non-linear part of this equivalent system; then a global resolution of the resulting polynomial system is achieved using a homotopy method. Our approach is compared to three existing methods: simple cell-mapping, a test exclusion based method and one relying on interval analysis. Section 2.1 presents a short overview and references for each existing method used for comparison while section 2.2 exposes the proposed method in a detailed way.

The four methods are applied to a Duffing oscillator and compared in section 3; this demonstrates efficiency of the suggested approach. Finally, its ability to treat more complex systems and the usefulness of such a global analysis algorithm is illustrated in section 4 through a simplified model of a rotor using a squeeze-film damper.

Global analysis methods at stake

This part is devoted to the four methods compared in the next section. First, the principles of three existing methods dedicated to global analysis are quickly explained. In a second part, the three main steps of the proposed method are detailed. Among the four treated methods, one (the cell-mapping) treats directly the differential system of equations (1) while the three others require its transformation into a non-linear (non-differential) system of equations [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF]. In the latter case the same transformation will be used for exclusion test methods, approaches relying on interval analysis and the proposed method. Finding all the T-periodic solutions of [START_REF] Alefeld | Interval analysis: theory and applications[END_REF] or finding all the solutions of ( 2) is equivalent if a proper transformation is used. A way to do so is exposed in part 2.2.1.

Ordinary differential equations.

M q + C q + Kq + f (q, q) = f e (t) [START_REF] Alefeld | Interval analysis: theory and applications[END_REF] where M, C and K are respectively mass, damping and stiffness matrices, f (q, q) stands for non-linear efforts depending on displacements and velocities and f e (t) is the excitation vector.

Non-linear (non-differential) system.

H x + Ĥ( x) = H e (2)
with x containing generalized unknowns, H denoting the linear part of equation, Ĥ( x) its non-linear part and H e its constant term.

Existing methods overview

First of all is the cell-mapping [START_REF] Hsu | Cell-to-Cell Mapping, a method of global analysis for non linear systems[END_REF] which consists in discretization of a bounded part of the state space in small "cells" (the rest of the state space is represented by a special "sink cell") and a short time integration of one or more points of each of these cells. The entire state space evolution along time can then be built and all kinds of solutions and their basins of attraction can be found. For example, let's consider a non-linear dynamical system of dimension 1 (with dimension 2 state space) and let q and q denote the displacement and velocity. Figure 1 shows a possible discretization used for the state space portion of interest (q, q) ∈ [q min , q max ] × [ qmin , qmax ]. The eight cells centers are startpoints for a time integration with duration τ (for non-autonomous system, τ often refers to excitation period, and to a characteristic time otherwise). The cell which contains the point after time-integration is the image cell of the cell in which time integration started. In figure 1, this mapping is drawn using arrows. On such a simple example, the analysis would lead to the detection of two cycles: one cycling through cells number 2,7 and 6 and the other one alternating between cells 4 and 5. Basins of attraction of these cycles are respectively cell 1 and cell 3. Cell 8 maps to the sink cell.

[Figure 1 about here.]

The illustrated branch of cell-mapping is the simple cell-mapping; some more elaborated versions can be found in literature ( [START_REF] Guder | Control of an adaptive refinement technique of generalized cell mapping by system dynamics[END_REF][START_REF] Ding | Application of an improved cell mapping method to bilinear stiffness aeroelastic systems[END_REF]). This method is quite natural because it deals with dynamical equations and seems attractive by finding all kinds of solutions and their basins of attraction, about what no other method can boast. The problem holds in the huge storage resource and computational time required to apply such a method.

Methods using exclusion tests require the transformation of ( 1) into [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF]. Then, a bounded part Ω of R ñ (where x varies) needs to be defined: that is where the solutions of (2) are sought. The algorithm tries to exclude portions of Ω using a cell-discarding-condition (CDC) [START_REF] Xu | A cell exclusion algorithm for determining all the solutions of a nonlinear system of equations[END_REF]: if the test is positive the cell contains no solution and is put away;

otherwise, the cell is kept and divided in two or more sub-cells that will be tested in turn. Figure 2 shows the general scheme of such methods. Their efficiency depends on both the quality of the test regarding the set of equations under study and its complexity: if it does not manage to eliminate big cells (big parts of Ω) it will fail in investigating the whole space in a reasonable amount of time, so it has to be of good quality but it should not cost too much (i.e. be too much complex) to avoid spending a lot of time when applied to a cell.

[Figure 2 about here.]

The authors chose to use a test based on phase I of the simplex method [START_REF] Yamamura | Interval solution of nonlinear equations using linear programming[END_REF]. This test is easy to implement and does not require any property for the non-linearity Ĥ( x). More specific tests can be found in [START_REF] Georg | Improving the efficiency of exclusion algorithms[END_REF][START_REF] Georg | A new exclusion test[END_REF].

Interval analysis based methods also process systems of form (2) rather than differential systems. As the previous method, they need a bounded part Ω of R ñ to be defined. This first cell is then partitioned leading to a collection of cells Γ. The main distinctive feature of these methods is that they no longer deal with real numbers but with intervals: usual operations such as "+", "-", etc. are redefined [START_REF] Alefeld | Interval analysis: theory and applications[END_REF]. Such algorithms try both to exclude cells and to reduce them (using other ways than simple splitting). Hansen and Walster [START_REF] Hansen | Global Optimization Using Interval Analysis[END_REF] propose an interesting algorithm that has two different approaches depending on the current cell size: some routines are dedicated to "big" cells and are cheap in term of computation cost; when a cell is "small" enough and has some properties, some other techniques, more costly are implemented to refine or exclude it. Figure 3 gives a simplified overview of this algorithm.

[Figure 3 about here.]

Main steps are Global approach processing (designed for big cells) and Local approach processing (designed for small cells). As it is a far complex algorithm, we do not give all the details for each main step but we rather give one representative routine for each. Let us denote X the cell being processed (interval of R ñ) and x a point of this cell (in R ñ).

• Global approach example: the hull consistency test. This method tries to take advantage of the linear part of system (2) which is the easiest part to deal with. Looking for a solution of (2), one can rewrite the equation in the following way ∀i ∈ {1, . . . , ñ}, H i j X j + Ĥi ( X) -H e i = 0

⇔ ∀i ∈ {1, . . . ,
H ik X′ k = -Ĥi ( X) + H e i ⊖ j k H i j X j , 1 ≤ i ≤ ñ with ⊖ standing for X ⊖ Y = [inf X -inf Y , sup X -sup Y].
Reusing already reduced Xk for reduction of Xm , k < m, this technique can lead to a cell exclusion (if Xk ∩ X′ k is empty) or at least to a great reduction using only simple operations.

• Local approach example: the newton like method. When the cell X is small enough, some expansion around a well chosen point x0 ∈ R ñ is a good approximation of non-linear part Ĥ( x) on the cell. Let's define

H( x) = H x + Ĥ( x) -H e . If H is regular enough ∀ x ∈ X, ∃ζ ∈ X, H( x) = H( x0 ) + J(ζ)( x -x0 ) with J( x) = ∂H i ∂ x j ( x) i j
.

The trick is to replace the unknown ζ by the whole cell X. Then if a solution x * exists in the cell X, on can write

0 ∈ H(x 0 ) + J( X)( x * -x0 )
This provides a new X′ interval containing all the solutions x of

J( X)( x -x0 ) = -H(x 0 )
which is a Newton like equation. As in the hull consistency test, X ∩ X′ may be empty or lead to X reduction.

Unfortunately, these three global analysis methods appear to be inefficient as soon as the system becomes greater than 2 or 3 degrees of freedom (dofs): if their mathematical background ensures to obtain all the solutions, the numerical resource and the computational time needed to treat the problem explode and most of the time one can predict neither this cost nor this time (see section 3).

Proposed polynomial homotopic method

A polynomial approximation of the non-linear part of system (2) is used in order to save time. Being in a welldefined space of multi-variable polynomial systems helps to apply more efficient methods like homotopies in order to find all the solutions of the algebraic system. In order to exhibit all T -periodic solutions of (1), one need to convert it into a non-differential equivalent system (2), then to turn this system into an algebraic one (i.e. a multi-variable polynomial system) and finally to solve this system "globally". Each of these three steps is detailed in the followings subsections.

From differential to non-differential equations

First of all, one needs to convert the dynamical formulation (1) into a non-differential system of form (2). This is done by making assumptions on the nature of the sought solution q(t) which, in the present case, is expected to be a T -periodic solution. One insists on the fact that such a step is usually followed by a local resolution algorithm (of Newton-Raphson type) which leads to one solution at most while a global resolution method will here be applied to the new set of equations in order to exhibit all its solutions.

The transformation relies on the harmonic balance method (HBM) [START_REF] Nacivet | A dynamic lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF][START_REF] Sinou | Center manifold and multivariable approximants applied to non-linear stability analysis[END_REF][START_REF] Nayfeh | Applied Nonlinear Dynamics[END_REF] which provides non-differential equations using a Galerkin approach. One assumes that q(t) is well enough described by a truncated Fourier series using N h harmonics:

q(t) ≈ a 0 √ 2 + N h k=1 a k cos(kωt) + b k sin(kωt)
with ω = 2π/T . It is straightforward to express q and q as a 0 , a k and b k functions and then to put back their expressions into (1), obtaining a set of n non-differential but time dependent equations with ñ = n(2N h + 1) unknowns (a 0 , a k , b k , 1 ≤ k ≤ N h ) and n denoting q size. The Galerkin approach generates the expected ñ equations by projecting this set of equations on each time dependent function used to decompose q(t) (i.e. 1/ √ 2, cos(kωt) and sin(kωt) with

1 ≤ k ≤ N h ).
The scalar product used is:

< f, g > T = 2 T T 0 f (t)g(t)dt
The resulting set of equations is obviously of form (2) (details are given in appendix A).

It is interesting to note that constant solutions are a special case of T -periodic solutions when q(t) is sought as a 0 only. Even if one mostly looks for cycles, it may be useful for some applications to be able to find all the fixed points solutions.

Reduction of the problem. Used in local search context, one often apply a Newton-Raphson scheme to these equations which means that a linear system has to be solved repeatedly: the smaller it is, the faster the solution is found. This assumption is all the more true in your case where computation cost grows severely with the system size; reducing the system can help saving hours for the resolution step. [START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF] proposes the exact condensation procedure that follows:

let us partition the unknowns xi which participate to the non-linear force expression Ĥ and those which do not and denote xnl the former and xl the latter. Rearranging (2), one gets:

          H ll H lnl H nll H nlnl                    xl xnl          +          Ĥl ( xnl ) Ĥnl ( xnl )          =          H el H enl         
The first rows, related to linear dofs, allow to write xl as a function of xnl :

xl = H -1 ll H el -Ĥ( xnl ) -H nl xnl
Using this new expression, the second row set becomes a non-linear equation in xnl of a much smaller size than ( 2):

H nlnl -H nll H -1 ll H nl xnl + Ĥnl ( xnl ) -H nll H -1 ll Ĥ( xnl ) = H enl -H nll H -1 ll H el ( 3
)
This equation is of size n(2N h + 1) where n denotes the non-linear dofs number. In many cases the structure is mostly linear and contains only a few non-linear elements. This makes use of this condensation interesting in local search as well as in global search context.

From a non-linear set of equations to a polynomial one

Non-linear system (2) or ( 3) has now to be "converted" into a polynomial system (or algebraic system) in order to apply the final global analysis method which relies on polynomial form of the system. The method will be exposed for the non condensed formulation (2).

Many ways of doing such an approximation exist. The authors chose to use a least squares method. This simple method may not be the most adequate but is easy to implement. A comparison with different methods should be achieved, dealing not only with the error resulting from the approximation but with the shape of the polynomials used and the consequences mainly on computational cost. Whatever the method may be it has to be noticed that only the non-linear part Ĥ( x) has to be fitted (by P): linear and constant parts have trivial polynomial expressions. Let us use the following notations : a) Polynomial system P and its components Pi :

P( x) =                  P1 ( x) . . . Pñ ( x)                  b)
Polynomial component Pi as a sum of J i coefficients and monomials products:

∀ i ∈ {1, . . . , ñ}, Pi ( xi ) = 1≤k≤J i β ik xα ik , β ik ∈ R, α ik ∈ N ñ, xα = ñ j=1 xα j j c
) Total degree d i of each component:

d i = deg( Pi ) = max 1≤k≤J i |α ik |, |α ik | = ñ j=1 α ik j
Using a least square method, one has to generate s points y ( j) = Ĥ( x( j) ) and to choose the monomials xα ik that will compose each component Pi of P. Monomials choice should be guided by the higher degree one need or can afford (see next section) to describe the true function Ĥ and obvious properties such as parity. Once this decision is made one finds minima of the cost functions φ i whose variables are coefficients β ik , 1 ≤ k ≤ J i (see Eq. ( 4)) by solving the linear system obtained when equating their derivatives to zero (Eq. ( 5) and ( 6)):

∀ i ∈ {1, . . . , ñ}, φ i (β i1 , . . . , β iJ i ) = s j=1 Pi ( x( j) ) -y ( j) 2 (4) ∀k ∈ {1, . . . , J i }, ∂φ i ∂β ik = 0 (5)
which is equivalent to:

A {β ir } 1≤r≤J i = b A rs = s j=1 ( x( j) ) α ir +α is , 1 ≤ r, s ≤ J i b r = s j=1 ( x( j) ) α ir y ( j) , 1 ≤ r ≤ J i ( 6
)
Possessing P, P which denotes the global polynomial system is obtained by simply adding the linear and constant part contributions:

P( x) =                  P 1 ( x) . . . P ñ( x)                  , P i ( x) = Pi ( x) + ñ j=1 H i j x j -H ei , 1 ≤ i ≤ ñ (7)
One should notice that using the polynomial approximation is often cheaper than using the original non-linear expression obtained by HBM: HBM method requires an AFT (Alternate Frequency-Time) procedure [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] for Ĥ term computation which is often costly enough. Using polynomials, one has to compute it s times only.

We will now focus on finding all zeros of this multivariate polynomial system.

Global resolution of the polynomial system using homotopies

Finding all the roots of a multivariate polynomial (in C) is not an easy task and one finds active research since 1977 with Garcia, Zangwill and Drexler work. More recently many papers dealt with solving multivariate polynomial systems using homotopy methods (see [START_REF] Wright | Finding all solutions to a system of polynomial equations[END_REF][START_REF] Li | The random product homotopy and deficient polynomial systems[END_REF][START_REF] Morgan | A homotopy for solving general polynomial systems that respects m-homogeneous structures[END_REF][START_REF] Morgan | Computing all solutions to polynomial systems using homotopy continuation[END_REF][START_REF] Zulehner | A simple homotopy method for determining all isolated solutions to polynomial systems[END_REF][START_REF] Watson | Probabiblity-one homotopies in computational science[END_REF] and [START_REF] Li | Numerical solution of multivariate polynomial systems by homotopy continuation methods[END_REF] for a general overview); the reference book [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF] written by Sommese and Wampler provides the full theory and great reflection about its numerical implementation.

Basically the homotopy method consists in a continuation from roots of a polynomial system Q to the roots of the system at stake P using the roots of the intermediate problem [START_REF] Hahn | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF] with λ varying from 1 to 0:

R( x, λ) = λQ( x) + (1 -λ)P, λ ∈ [0, 1] (8)
Having all complex zeros of P, it is easy to extract the real ones which are approximate solutions of the initial nonlinear system (2).

Obviously Q should have "easy-to-find" roots and it should also match some conditions to ensure that the homotopy process will provide all of P zeros. First of all it should have as many zeros as P because each P zero must be linked to a Q zero through a continuous path. Then, it must share structural properties with P so as to avoid bifurcation along continuation (which would make the whole procedure inefficient). A simple Q that respects these conditions and has easy to find roots is:

Q( x) = γ                  xd 1 1 -1 . . . xd ñ ñ -1                  , γ ∈ C, d i = deg(P i ) (9)
A total-degree-homotopy theorem ([21, Chap. 8, p. 123]) ensures that following all Q roots, one will find all the nonsingular solutions of P( x) = 0 (solved in C ñ). The problem is that most of the N = ñ i=1 d i paths followed will not converge to a sought solution but diverge to infinity and so, are worthless to follow. If dealing with divergence is no real problem (solved using a homogeneization technique explained below), following a huge number of paths can make this method become unusable practically. That is why most of the papers deal with building initial polynomial sets Q that have the fewest roots while still having easy to find ones. Among these methods are the total degree homotopy method and the multihomogeneous homotopy method (often referred to as m-homogeneous homotopy). The first one consists in using a polynomial Q as expressed in [START_REF] Hansen | Global Optimization Using Interval Analysis[END_REF] considering that all variables xi play the same role: no real analysis of P is required but it generates the maximum number of paths (majored by Bézout number N). The second one uses m groups of variables taking into account, for example, that most of them only have a linear contribution while a few have a non-linear (high degree) contribution.

The number of paths generated using such a method is given by the Bézout theorem adapted to multi-homogeneous structures [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF]; optimal partitioning methods can be found in [START_REF] Wampler | Bézout number calculations for multi-homogeneous polynomial systems[END_REF][START_REF] Liang | Computing the optimal partition of variables in multi-homogeneous homotopy methods[END_REF][START_REF] Tao | A hybrid algorithm for multi-homogeneous bzout number[END_REF]. Generalizing this, one can use more than one homogeneization variable. For example m variables denoted z1 , . . . , zm .

If generalized unknowns are partitioned in m groups

Z k = { xi , i ∈ K k }, with ∪ 1≤k≤m K k = {1, . . . , ñ} and ∀1 ≤ i < j ≤ m, K i ∩ K j = ∅,
a homogeneization variable zk can be affected to each group of generalized unknowns Z k .

A degree table D = (d i j ) 1≤i≤ñ,1≤ j≤m can then be established with d i j denoting the maximum degree of group K j variables observed in monomials of P i . Then, for each P i component of P a multi-homogeneization step can be achieved by bringing each P i monomial to the same degree j d i j using the m homogeneization variables: each generalized unknown xi in group Z k is temporarily replaced by xi /z k and the resulting denominators are cleared by multiplying P i by 1≤k≤m zd ik k . The multi-homogenized polynomial P h (z 1 , . . . , zm , x1 , . . . , xñ ) has then component P hi whose monomials share the same degree j d i j but where generalized unknowns of different groups are treated differently. This approach is interesting whenever some variables occur with low degrees only (generalized unknowns related to linear dofs) and a few of them occur with high degrees (generalized unknowns related to non-linear dofs).

an initial polynomial Q as specified in [START_REF] Hsu | Cell-to-Cell Mapping, a method of global analysis for non linear systems[END_REF] (the important point is that P and Q have the same degree table D) can then be used and has easy to find roots.

Q( x) = γ                  ( xd 11 K1 -1) × • • • × ( xd 1m Km -1) . . . ( xd ñ1 K1 -1) × • • • × ( xd ñm Km -1)                  , γ ∈ C (10)
where xKi stands for one of the generalized unknowns in group K i .

These homotopy methods have a great advantage over cell-exclusion or interval analysis based ones: they allow computational time estimation. Knowing in advance the number of paths to be followed, it is easy to compute a few of them, obtaining then a rough estimation for the whole stack. Plus, they do not require to be given an arbitrary initial bounded space where to look for the solutions whereas the three others do: they ensure to get all the solutions, letting us decide a posteriori if they have a physical meaning or interest, or not.

Comparison with other global analysis methods on a simple Duffing oscillator

In order to evaluate the efficiency of the proposed approach, it will be compared to three different methods using a classical non-linear Duffing oscillator. The results exposed come from the PhD Thesis [START_REF] Sarrouy | Analyse globale de systmes mcaniques non-linaires -Application la dynamique des rotors[END_REF] where details can be found.

The forced Duffing oscillator at stake is : q + 2ξω 0 q + ω 2 0 q + βq 3 = Γ cos(ωt) [START_REF] Li | Numerical solution of multivariate polynomial systems by homotopy continuation methods[END_REF] with numerical values ξ = 0.02, ω 0 = 1 rad.s -1 , β = 1 s -2 and Γ = 0.1 m.s -2 . This forced oscillator has the frequency response given in figure 4; this curve is obtained using a continuation method. While looking for multiple solution sets, global analysis methods are applied for ω = 1.2 rad.s -1 : three solutions -two stable S 1 and S 3 and one unstable S 2 -coexist as shown on the figure.

[Figure 4 about here.]

The main characteristics of used methods are given below:

(a) simple cell-mapping: bounded state space portion (q, q) considered is [-2, 2] × [-3, 3] using 100 divisions in each directions (leading to 10 000 cells to be examined). Each cell center is a startpoint for a time-integration scheme with duration τ = 2π/ω.

(b) exclusion test method: dynamic equation ( 11) is turned into a non-linear system using HBM with different harmonic numbers N h in order to obtain different sizes for the problem (2). The exclusion test used is based on simplex method [START_REF] Yamamura | Interval solution of nonlinear equations using linear programming[END_REF]. (d) polynomial homotopy method: dynamic equation ( 11) is turned into a non-linear set using HBM with different harmonic numbers N h in order to obtain different sizes for the problem ( 2). This problem is already polynomial (total degree 3): no polynomial approximation step is performed. A total degree homotopy method is used to build the initial polynomial Q.

Each
Every method is implemented by using Matlab 7; the calculations are made on a Core 2 Duo E6600 (4Go RAM).

[Figure 5 about here.]

First, methods are applied to the Duffing oscillator, with various harmonic decompositions for methods (b) to (d) in order to get some reference times. Then, extrapolations are computed to investigate the methods efficiency for larger systems. Figure 5 summarizes the results. As cell-mapping is applied directly on the dynamic equation while other methods are applied to the transformed non-linear set of equations, results are not exposed the same way. The first graphic deals with the simple cell-mapping method. It shows the computation time required to process systems of 1 to 3 dofs using 10, 20, 30, 40, 50 or 100 divisions in each direction of the investigated bounded portion of the state space. The second graphic gives the time required for systems with up to 10 xi unknowns (with a degree 3 polynomial approximation and an initial polynomial based on total degree in case of method (d)). The number of dofs related to this number of generalized unknowns depends on the harmonic decomposition used: 5 dofs using only 1 harmonic (no constant) leads to 10 generalized unknowns as do 2 dofs with a constant term and 2 harmonics. Axis of ordinate is logarithmic and representative time units are indicated. This shows that the proposed method is obviously faster than methods (b) and (c); for comparison with method (a), one can consider a 2 dofs system with a constant term and 2 harmonics decomposition (10 xi unknowns): the cell-mapping with 50 divisions in each direction and the polynomial homotopy method seem then cost equivalent. Drawbacks of the cell-mapping method is the need of a bounded state space where solutions are sought and the arbitrary number of divisions used: with only 50 divisions, one may be brought to refine some areas of interest and so to increase the time required to get the solutions. Furthermore, required time growth with system size is much bigger for the cell-mapping than for the polynomial homotopy method. Finally, the proposed approach can be easily parallelized as each path can be followed independently.

According to its efficiency and potential, one proposes to evaluate this most promising technique on a some more complex example.

Application to a squeeze-film supported rotor

We will now apply homotopy methods to a simple model of rotor damped using a squeeze-film bearing (fig. 6).

[Figure 6 about here.]

This system has already been studied in [START_REF] Zhao | Eccentric operation and blade-loss simulation of a rigid rotor supported by an improved squeeze film damper[END_REF] using local methods. It consists in a stiff and damped rotor m 1 , c 1 , k 1 (including rolling bearing and squeeze-film inner ring) mounted on a squeeze-film bearing with outer ring m 2 linked to the "basement" through c 2 , k 2 ; excitation comes from the unbalance U. Dynamic equations are (12) (using notations of appendix B) where the non-linearity is induced by the squeeze-film through F x and F y efforts:

M q + C q + Kq +                          -F x -F y F x F y                          =                          U cos τ U sin τ 0 0                          (12)

From dynamic to non-linear equations

Using the classical centered circular synchronous orbiting assumption for the rotor and a concentric movement assumption for the whole system, one can simplify the equations by rewriting q and the non-linear efforts F x and F y :

q(t) =                          a 1 b 1 a 2 b 2                          cos τ +                          b 1 -a 1 b 2 -a 2                          sin τ (13)
F r = - B s 2 ǫ 2 (1 -ǫ 2 ) 2 ω ( 14
)
F t = - B s π 2 ǫ (1 -ǫ 2 ) 3/2 ω (15)
As harmonic coefficients in cos τ and sin τ are bound because of circular orbiting assumption, one only needs to solve equations obtained by balancing [START_REF] Li | The random product homotopy and deficient polynomial systems[END_REF] with respect to cos τ:

                    K -M + C                     0-1 1 0 0-1 1 0                                         H x +                          -H x -H y H x H y                          Ĥ( x) =                          U 0 0 0                          H e , H x/y ( x) = 1 π 2π 0 F x/y ( x) cos(τ)dτ (16)

Building the polynomial approximation

We now need to build a polynomial approximation P( x) of the non-linear part Ĥ( x). First let us notice that it depends only on ǫ that is the relative displacements X 2 -X 1 and Y 2 -Y 1 : the change of variables ( 17) is applied.

q = Bq r , B =                          1 1 1 1 1 1                          (17)
From now on x will denote harmonics vector of q r decomposition:

q r =                          a 1 b 1 a r b r                          x cos τ +                          b 1 -a 1 b r -a r                          sin τ (18)
This change of variables has small changes on ( 16) but allows us to use a 2 variables polynomial (whose variables are a r and b r ) instead of a 4 variables one which saves monomials used to express the polynomial and so computation time required to evaluate it and even more, as we will see in the next subsection when applying a homotopy method.

[Figure 7 about here.]

As mentioned in section 2.2.2, a least squares method is used to fit Ĥ. Two polynomial approximations are built: one using degree 6 polynomials and the other using degree 12 polynomials. The choice is made to generate points y ( j) used for fitting that respect 0 ≤ ǫ ≤ 0.6; that is, a r and b r couples are picked up such that 0 ≤ a 2 r + b 2 r ≤ 0.6. This is motivated by the physical reality: formulas used to approximate squeeze-film efforts are not valid beyond such eccentricities.

Figure 7 allows some qualitative comparison between the true values of H x and the ones obtained using polynomial approximations while figure 8 gives quantitative results for degree 12 approximation of H x : mean relative error is plotted for small rings along ǫ and shows good agreement with the original function with less than 5% error for almost all ǫ between 0.2 and 0.6. Bad results obtained for very small ǫ values can be explained by almost zero value of H x for small eccentricities -such eccentricities are not of great interest anyway because non-linear effects are then almost negligible.

[Figure 8 about here.] Polynomial systems obtained after fitting are augmented with linear and excitation parts H and H e to obtain P( x) as defined in sec. 2.2.2.

Global analysis using homotopies

Different homotopy methods were tested for one configuration (s = 1) and for each polynomial approximation (c) Total degree homotopy on reduced system: after applying the exact condensation of section 2.2.1, a new polynomial system is computed (containing the only 2 non-linear variables a r and b r ) and a total degree homotopy method is applied, generating d 2 paths.

The results in terms of number of paths and time required to obtain the complex solutions are summed up in table 1.

[Table 1 about here.] Every method ends in a reasonable amount of time but there exists a great gap between each of them, the best being the total degree homotopy used to solve the condensed system. This method was then used with degree 12 approximation for several s values. The real solutions -the only ones of interest -are kept and plotted using red points on figure 9.

[Figure 9 about here.]

The solid blue and dashed green curves were obtained using a continuation scheme (see [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF] for details). The solid (lower) one is the classical dynamic response obtained when using a continuation scheme starting from a local search for a small s value. The dashed (upper) one is the curve obtained after a local search using one of the global analysis results as a startpoint. The gap between the second branch depicted by global analysis points and the one obtained using a local continuation scheme results from the error made by the polynomial approximation for large eccentricities.

This shows how important it is to possess a global analysis tool: without it, one would not look for an another branch of solutions, disconnected from the solid one, appearing without any obvious bifurcation.

A stability study is conducted for these two branches using Floquet theory (figure 10): the second branch has stable solutions of large eccentricities coexisting with solutions of the classical one. Figure 11 shows that the second branch, on its stable part, consists in orbits of quite the same magnitude than orbits of the first branch for the rotor, but squeeze-film orbits are much greater; in fact, the squeeze-film outer-ring orbits outside the journal orbit.

[Figure 10 One already mentioned that the π-film model was not a good one to describe such large eccentricities so squeezefilm damper users should not be afraid of this new branch until other tests using proper models are done. Nevertheless, this points out the need to be careful when working with such great eccentricities, and shows the ability of the method to exhibit unexpected disconnected branches of solutions, which was at stake.

Conclusion

This paper exposes a global analysis method able to find all the T -periodic solutions of a dynamical non-linear system. It is based on three major steps: first a harmonic balance to transform dynamic equations into a non-linear non-differential system; this system is then put into a polynomial form in order to proceed to the global analysis (finding all its roots) using a homotopy technique. It has be shown to be the most efficient approach when compared to three other global analysis methods applied to a Duffing oscillator.

A simple model of rotor supported by a squeeze-film damper emphasizes then the ability of the method to find disconnected solutions and its computational efficiency. High degree polynomial approximations of non-linear part can be treated in a reasonable amount of time which allows to treat complex non-linearity expressions.

The harmonic balance was applied using a known period T ; further work will be achieved so as to let it free, by adding for example a phase equation [START_REF] Seydel | From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis[END_REF]. Another interesting improvement will be to compare different polynomial approximation methods linked with variables space partitioning to save computational time. Finally, parallelization could be easily implemented as each path of the homotopy step can be processed independently which would help to save a lot of computation time.

A. HBM elements definition a) Vector of generalized unknowns

x: x = a 0 , a 1 , b 1 , . . . , a N h , b N h t ∈ R n(2N h +1) Displacement q(t) is a function of x: q(t) = φ ω ( x) = a 0 √ 2 + N h k=1 a k cos(kωt) + b k sin(kωt) b) Linear part H x: H =                           Λ 0 Λ 1 . . . Λ N h                           , with                  Λ 0 = K ∈ M N (R) Λ k =           K -(kω) 2 M -(kω)C (kω)C K -(kω) 2 M           , Λ k ∈ M 2N (R), 1 ≤ k ≤ N h c) Non-linear part Ĥ( x) and constant term H e : Ĥ( x) =                                    < f (φ ω ( x)), 1/ √ 2 > T . . . < f (φ ω ( x)), cos(kωt) > T < f (φ ω ( x)), sin(kωt) > T . . .                                    1≤k≤N h , H e =                                    < f e (t), 1/ √ 2 > T . . . < f e (t), cos(kωt) > T < f e (t), sin(kωt) > T . . .                                    1≤k≤N h B. Squeeze-film example notations B.1. Initial system of equations                          m 1 x ′′ 1 +c 1 x ′ 1 +k 1 x 1 = f x +m 1 uω 2 cos(ωt) m 1 y ′′ 1 +c 1 y ′ 1 +k 1 y 1 = f y +m 1 uω 2 sin(ωt) m 2 x ′′ 2 +c 2 x ′ 2 +k 2 x 2 =-f x m 2 y ′′ 2 +c 2 y ′ 2 +k 2 y 2 =-f y
, with z ′ = ∂z ∂t , leads to [START_REF] Li | The random product homotopy and deficient polynomial systems[END_REF] after division by m 1 Cω 2 . [START_REF] Ding | Application of an improved cell mapping method to bilinear stiffness aeroelastic systems[END_REF]. SF contribution using a short bearing approximation and a π-film cavitation model

B.2. Dimensionless variables

τ = ωt, ż = ∂z/∂τ, s = ω/ω 1 , α = m 2 /m 1 , β = k 2 /k 1 , ξ 1 = c 1 /(m 1 ω 1 ), ξ 2 = c 2 /(m 2 ω 2 ), X i = x i /C Y i = y i /C, U = u/C, F x = f x /(m 1 Cω 2 ), q = {X 1 , Y 1 , X 2 , Y 2 } t , ǫ = (X 1 -X 2 ) 2 + (Y 1 -Y 2 ) 2 B.3. Matrices M = diag(1, 1, α, α), C = 1 s diag(ξ 1 , ξ 1 , ξ 2 √ αβ, ξ 2 √ αβ), K = 1 s 2 diag(1, 1, β, β) B.
F r = - B s π 2 1 + 2ǫ 2 (1 -ǫ 2 ) 5/2 ǫ + 2 ǫ 2 (1 -ǫ 2 ) 2 φ , F t = - B s 2 ǫ (1 -ǫ 2 ) 2 ǫ + π 2 ǫ (1 -ǫ 2 ) 3/2 φ , F x = F t cos φ + F r sin φ, F y = -F r cos φ + F t

  To better understand the difference between these two methods, one needs to describe what homogeneization means. The shortest way to explain it is: turning each monomial with degree d i of a polynomial with total degree d = max d i into monomials with the same degree d by multiplying them by x(d-d i ) 0 where x0 is a homogeneization variable. Then ∀( x0 , x1 , . . . , xñ ) ∈ C (P+1) , ∀µ 0, P(µ( x0 , x1 , . . . , xñ )) = µ d P( x0 , x1 , . . . , xñ ) which is at least an interesting property for dealing with divergence: if P( x1 , . . . , xñ ) is the polynomial system to solve, let's denote P h ( x0 , x1 , . . . , xñ ) the equivalent system where each P k component is homogeneized. Adding an unknown x0 , one needs to add an equation; for example, max 0≤i≤P | xi | = 1 This leads to a homogeneized continuation system R h (( x0 , . . . , xñ ), λ) = 0. If during continuation, one of the xi become huge to satisfy R h (( x0 , . . . , xñ ), λ) = 0, one can use µ = 1/| xi | and rewrite: R h (µ( x0 , . . . , xñ ), λ) = µ d R h (( x0 , . . . , xñ ), λ) = 0 µ( x0 , . . . , xñ ) is not that huge anymore and still satisfies the continuation equation. Plus, solutions obtained when λ tends towards 0 whose x0 component tends towards 0 are points to infinity.

  degree d ( d = 6 or d = 12): (a) Total degree homotopy on full system: the ñ = 4 generalized unknowns are kept and considered the same way. Each of the four equations P i ( x) is of degree 6 or 12, leading to d 4 paths to follow. (b) 2-homogeneous homotopy on full system: the 4 unknowns are partitioned into 2 subsets; one of them contains linear variables a 1 and b 1 and the other contains non-linear ones a r and b r . The number of paths to continue is then 6d 2 .

  about here.] [Figure 11 about here.]

sin φ B. 5 .Figure 1 : 1 Figure 2 :Figure 3 :Figure 4 :Figure 5 :k 1 , c 1 k 1 , c 1 k 2 , c 2 k 2 2 φFigure 6 : 1 - 1 L 2 Figure 7 :Figure 8 :Figure 9 :Figure 10 : 1 -Figure 11 :

 511234511112222611278910111 Figure 1: Simple cell-mapping example

  ñ}, ∀k ∈ {1, . . . , ñ}, H ik Xk +

	j k	H i j X j + Ĥi ( X) -H e i = 0
	with ñ denoting x dimension. It is then obvious that each component interval Xk can be reduced to its intersection
	with X′ k satisfying

List of Figures


E Sarrouy

A Grolet

F Thouverez
email: fabrice.thouverez@ec-lyon.fr

Global and bifurcation analysis of a structure with cyclic symmetry

Keywords: Global analysis, Non-linear, Homotopy, Rotordynamics, Cyclic symmetry

This article combines the application of a global analysis approach and the more classical continuation, bifurcation and stability analysis approach of a cyclic symmetric system. A solid disc with four blades, linearly coupled, but with an intrinsic non-linear cubic stiffness is at stake. Dynamic equations are turned into a set of non-linear algebraic equations using the Harmonic Balance Method. Then periodic solutions are sought using a recursive application of a global analysis method for various pulsation values. This exhibits disconnected branches in both the free undamped case (non-linear normal modes, NNMs) and in a forced case which shows the link between NNMs and forced response. For each case, a full bifurcation diagram is provided and commented using tools devoted to continuation, bifurcation and stability analysis.

Introduction

In this paper, both free and forced vibrations of a non-linear cyclic symmetric structure are studied.

The structure is composed of n b identical substructures which undergo large strains. This system is typical when one studies bladed disks [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Vakakis | A multiple-scales analysis of nonlinear, localized modes in a cyclic periodic system[END_REF]. After modeling the system, a set of coupled non-linear differential equations in which non-linearity appears by cubic terms is obtained .

In the linear case, the study of the linear normal modes (LNMs) reveals a majority of double eigenfrequencies, corresponding to distinct eigen-forms [START_REF] Samaranayake | Subharmonic oscillations in harmonically excited mechanical systems with cyclic symmetry[END_REF]. These eigenforms are associated with nodal diameter vibration modes. As these LNMs arise from an eigenvalue problem, there are always as many LNMs as degrees of freedom (dofs). Moreover, in the free or forced case, only one solution exists for a given frequency.

In the case of non-linear systems with cyclic symmetry, it has been shown that the number of nonlinear normal modes (NNMs, [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]) can exceed the number of dofs, the extra NNMs being generated through bifurcations or internal resonances [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Vakakis | Nonlinear normal modes and their application in vibration theory: an overview[END_REF].

Localized non-linear modes are an example of this property; they correspond to a free motion in which only a few substructures vibrate with non negligible amplitude and they have no counterpart in the linear theory [START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Vakakis | A very complicated structure of resonances in a nonlinear system with cyclic symmetry: non linear forced localization[END_REF]. In the forced case, these additional NNMs give rise to a number of additional resonances leading to multiple solutions as in [START_REF] Vakakis | A very complicated structure of resonances in a nonlinear system with cyclic symmetry: non linear forced localization[END_REF] where Vakakis showed a very complicated structure of resonance for a non-linear system with cyclic symmetry by using the multiple scales method. He showed that for a given frequency many solutions can coexist, some of them being stable.

Not only can multiple solutions coexist but also can they be disconnected from each other. In this latter case, classical methods based on continuation and bifurcation analysis fail at finding the disconnected branches of solutions. If one wants to dimension a structure properly -by considering all the possible solutions -one then needs to adopt a global analysis (GA) approach. Several solutions are mentioned in the literature (see [START_REF] Sarrouy | Global search of nonlinear systems periodic solutions: A rotordynamics application[END_REF] for an overview) with a common drawback even for small systems which is the computation cost. The GA method proposed in this paper takes advantage of the cubic form of the non-linearity combined with a reformulation of equations through the harmonic balance method (HBM); the resulting systems in the free undamped case for searching NNMs as well as in a forced case can then be solved (globally) in a reasonable amount of time.

Section 2 describes the system and its dynamical and HBM equations. The global analysis principle is then explained in section 3.1; this is followed by recalls on continuation methods, bifurcation and stability analysis in section 3.2. Finally these methods are applied to the undamped free system in order to find NNMs in section 4 and to a forced case which exhibits a very rich response in section 5.

Studied system

General description and dynamical equations

The studied structure has a cyclic symmetry property and can therefore be broken up in n b identical sectors (Fig. 1). Each sector is modeled by a thin rectangular plate clamped to the rigid disk which itself is fixed (Fig. 2). Consecutive plates are coupled by a linear stiffness while nonlinearity is introduced by taking into account their large deflection. Plane stress assumption as well as the Love-Kirchhoff hypothesis (cross-sections exhibit solid body motion and remain perpendicular to the deformed surface of the middle sheet) are made. Then, plates displacements are entirely parameterized by their middle sheet transverse displacement w j , 1 ≤ j ≤ n b . Moreover the material is assumed to follow a standard bi-dimensional Hooke's law leading to the following expression for the strain energy of plate j:

U j = 1 2 Eh (1 -ν 2 ) Lx x=0 Ly/2 y=-Ly/2 1 2 ( ∂w j ∂x ) 2 2 + 1 2 ( ∂w j ∂y ) 2 2 + 2ν 1 2 ( ∂w j ∂x ) 2 1 2 ( ∂w j ∂y ) 2 + 1 -ν 2 ∂w j ∂x ∂w j ∂y 2 dxdy + 1 2
Eh 3 12(1 -ν 2 ) Lx x=0 Ly/2 y=-Ly/2 ∂ 2 w j ∂x 2 2 + ∂ 2 w j ∂y 2 2 + 2ν ∂ 2 w j ∂x 2 ∂ 2 w j ∂y 2 + 1 -ν 2 2 ∂ 2 w j ∂x∂y 2 dxdy ( 1
)
where E is the Young's modulus and ν is the Poisson's ratio.

The energy V j of the linear stiffness localized at (x r , y r ) = (L x /4, 0) between plates j and j + 1 is given by:

V j = 1 2 k(w j (x r , y r ) -w j+1 (x r , y r )) 2 (2) for 1 ≤ j ≤ n b with convention j + 1 = 1 if j = n b .
By neglecting rotary inertia the kinetic energy T j of plate j is given by :

T j = 1 2 ρh Lx x=0 Ly/2 y=-Ly/2 ẇj 2 dxdy (3)
In this paper, only an harmonic force, orthogonal to the plate, localized in (x f , y f ) = (L x , 0) is considered. The work W j due to such an excitation on plate j is given by:

W j = w j (L x , 0)F ej (4)
The total energies U , T , V and W are then given by the sum over the number of plates n b of the different local energies U j , T j , V j and W j . Equations of motions are finally derived by using Lagrange's equations along with a Rayleigh-Ritz approximation [START_REF] Liew | pb2-rayleigh-ritz method for general plate analysis[END_REF].

[Figure 1 x = 0) and leads to the following interpolation for transverse displacement of the j-th blade:

w j (x, y) = q j Φ(x, y) for 1 ≤ j ≤ n b (5)
By applying Lagrange's equations, and by adding a damping term, the following motion equations are obtained:

[M] Ẍ + [C] Ẋ + [K]X + βX 3 = F e (t) (6)
with the notations X = (q j ) 1≤j≤n , and X 3 = (q 3 j ) 1≤j≤n . The vector F e = F e0 cos(ωt) stands for the external forces amplitude, [M] = [I] is the mass matrix, [C] = δ[I] is the damping matrix, [K] is the stiffness matrix given by:

[K] =         α + 2c -c 0 -c -c α + 2c -c 0 0 -c α + 2c -c -c 0 -c α + 2c         (7)
and β is the non-linear stiffness coefficient. The definitions of α, β, c, and δ can be found in A along with their numerical values.

Harmonic balance method

The harmonic balance method (HBM) is widelyused for the study of non-linear systems. Numerous applications can be found in the literature, showing its ability to treat strongly non-linear systems like friction between blades and casing [START_REF] Laxalde | Qualitative analysis of forced responce of blisks with friction ring dampers[END_REF][START_REF] Laxalde | Complex nonlinear modal analysis for mechanical systems: application to turbomachinery bladings with friction interfaces[END_REF] or geometric non-linearities [START_REF] Lewandowski | Computational formulation for periodic vibration of geometrically nonlinear structurespart 1:theoretical background[END_REF][START_REF] Ribiero | Nonlinear vibration of plates by the hierarchical finite element and contination method[END_REF]. One major advantage of the method is that it requires no assumption about the non-linearities magnitudes and uses the same procedure for strongly and weakly non-linear models.

The HBM consists in a decomposition of the solution X in a truncated Fourier series:

X(t) = A 0 + N h k=1 A k cos(kωt) + B k sin(kωt) (8)
Injecting this development [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF] in equation [START_REF] Huber | Polyhedral end games for polynomial continuation[END_REF], and by projecting equations on the [1, (cos(kωt), sin(kωt)) 1≤k≤N h ] basis using the following scalar product

f, g = 2π ω 0 f (t)g(t)dt (9)
one gets a system of ñ = n × (2N h + 1) non-linear algebraic equations with n×(2N h +1)+1 unknowns

A k , B k and ω.
The number of harmonics retained N h is a very important parameter. Generally, the higher N h is, the better the solution. However, in the case where the number of harmonics selected is high, the solution procedure can quickly become difficult and time consuming. Fortunately, in most cases the series converges fast enough and leads to systems with reasonable dimensions. In this article only the first harmonic is going to be retained (N h = 1) and the constant term A 0 is dropped due to the symmetry of the system. Depending on whether the system is free or forced, different formulations can be obtained.

Forced case

The equation to be solved in the forced case is given by [START_REF] Huber | Polyhedral end games for polynomial continuation[END_REF]. The solution X is sought of the following form:

X(t) = A cos(ωt) + B sin(ωt) (10)
which leads to the following set of 2n algebraic equations:

[H l ] X + H nl ( X) = H e (11)
where

X =    A B   
is the unknown in the frequency domain; it follows that ñ = 4 × 2 = 8. [H l ] corresponds to the linear part of the HBM:

[H l ] =   [K] -ω 2 [M] ω[C] -ω[C] [K] -ω 2 [M]   (12)
H e corresponds to the constant part related to the excitation

H e =    F e0 0    (13)
and H nl corresponds to the non-linear part

H nli ( X) = 3 4 β(A 3 i + A i B 2 i ) H nln+i ( X) = 3 4 β(B 3 i + B i A 2 i ) for 1 ≤ i ≤ n (14)
Notice that in this case the solution frequency ω is taken equal to the one of the excitation force.

Free Case

The equation to be solved in the free undamped case is given by:

[M] Ẍ + [K]X + βX 3 = 0 (15)
and the solution X is sought of the following form:

X(t) = A cos(ωt) (16)
The fact that only the cosine terms are retained corresponds to a phase condition in which all initial velocities are set to zero. With this approximation for the solution X, ñ = 4 and the following set of algebraic equations is obtained:

([K] -ω 2 [M])A + 3 4 βA 3 = 0 ( 17
)
As there is no excitation force, the free solution frequency ω is also an unknown.

Theory

This section is devoted to theory on both global analysis (GA) aspects and continuation coupled to stability and bifurcation analysis.

Global analysis

When aiming at finding all the solutions of a nonlinear system, looking at the whole state space seems the right way to proceed. This is what most global analysis methods such as cell-mapping [START_REF] Hsu | Cell-to-Cell Mapping, a method of global analysis for non linear systems[END_REF],

cell exclusion algorithm [START_REF] Xu | A cell exclusion algorithm for determining all the solutions of a nonlinear system of equations[END_REF] or algorithms based on interval analysis [START_REF] Hansen | Global Optimization Using Interval Analysis[END_REF] do. The main drawback of all these methods is their huge cost with regards to their computational resource consumption (see [START_REF] Sarrouy | Global search of nonlinear systems periodic solutions: A rotordynamics application[END_REF] for a comparison). In order to bypass this problem, one can take advantage of the algebraic form of our problem formulated in the frequency space and use a dedicated method to solve it globally relying on homotopy [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF][START_REF] Sarrouy | Global search of nonlinear systems periodic solutions: A rotordynamics application[END_REF][START_REF] Sarrouy | Analyse globale de systèmes mécaniques non-linéaires -Application à la dynamique des rotors[END_REF][START_REF] Li | Numerical solution of multivariate polynomial systems by homotopy continuation methods[END_REF]. Indeed, the nonlinearity introduced being cubic, equations [START_REF] Laxalde | Qualitative analysis of forced responce of blisks with friction ring dampers[END_REF] or ( 17) can be considered as polynomials in X variable. We will now denote P the corresponding polynomial.

The main idea of this method is to embed the polynomial P to solve in a space whose parameters are the polynomial coefficients. Then, starting from a polynomial Q whose zeros (in C ñ) are easy to compute, a continuation scheme is applied to slide smoothly from these zeros to the ones of interest by substituting slowly P to Q which is, namely, applying a homotopy method. The system solved is

R( X, λ) = λQ( X) + (1 -λ)P( X), λ ∈ [0, 1] (18)
with λ varying from 1 to 0.

In order to be able to guarantee that all P zeros will be reached, one has to respect a few rules: first of all, as each non-singular zero of P must be linked to a non-singular zero of Q, the initial polynomial must have at least as many zeros as P. Then, the paths linking Q zeros to P zeros must be smooth (undergo no bifurcation when λ slides from 1 to 0) which implies that both polynomials must share structural properties [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF]Chap. 8]. A basic polynomial which has easy to find zeros and respects these rules is

Q( X) = γ          Xd1 1 -1 . . . Xdñ ñ -1          , γ ∈ C, d i = deg(P i ) (19)
Such a polynomial has N = ñ i=1 d i zeros, that is it has the maximum number of zeros a polynomial with such total degrees (maximum degree d i among its monomials for each component P i ) can have. A homotopy process with this kind of initial polynomial is called a total degree homotopy.

The total degree homotopy theorem demonstrated in [25, chap. 8, p. 123] ensures that the N zeros of Q will lead to all the non-singular zeros of P. In our frame, only zeros with real components will be considered.

In many cases, P has only a few complex zeros. Starting with a lot more means a waste of time following paths that will mostly diverge towards infinity (as a homogeneization technique is used, they will remain numerically bounded, see [START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF][START_REF] Wright | Finding all solutions to a system of polynomial equations[END_REF]). That is why, in most cases it is interesting to work on the starting polynomial to reduce the number of starting zeros while still respecting the rules enunciated previously. Most current work aims at exhibiting such systems and automatizing their construction ( [START_REF] Morgan | A homotopy for solving general polynomial systems that respects mhomogenous structures[END_REF][START_REF] Morgan | Coefficient-parameter polynomial continuation[END_REF][START_REF] Wampler | Bézout number calculations for multihomogeneous polynomial systems[END_REF][START_REF] Sommese | The Numerical Solution of Systems of Polynomials Arising in Engineering and Science[END_REF][START_REF] Huber | Polyhedral end games for polynomial continuation[END_REF][START_REF] Lee | Hom4ps-2.0: a software package for solving polynomial systems by the polyhedral homotopy continuation method[END_REF]). As will be shown later in this paper, such work is useless in the studied case. Using polynomial Q defined in [START_REF] Nayfey | Applied nonlinear dynamics[END_REF], 3 8 = 6561 starting zeros will be generated in the forced case (HBM equation [START_REF] Laxalde | Qualitative analysis of forced responce of blisks with friction ring dampers[END_REF]) and only 3 4 = 81 paths will be followed in the free one (HBM equation ( 17)).

Continuation, bifurcation and stability

In this section, the classical tools to perform continuation and bifurcation analysis of fixed points (solutions in frequency space are constant) as well as the stability analysis of periodic solutions relying on Floquet coefficients evaluation are briefly recalled.

Continuation algorithm

A general system of equations used by a continuation algorithm is

G(Y, µ) = 0 ( 20
)
where G is a non-linear function, Y is the vector of unknowns and µ is the continuation parameter. In the present case Y = X, the HBM unknowns and µ = ω, the free or forced frequency. The continuation procedure builds the path solution of [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF] for some µ range by iterating two steps. First, a point is predicted using one or more of the previous points of the path, yielding an estimation of the next point. This predicted solution is then corrected iteratively until convergence is achieved.

In this study an arc-length continuation is used:

G is obtained by augmenting the HBM set of equations [START_REF] Laxalde | Qualitative analysis of forced responce of blisks with friction ring dampers[END_REF] or [START_REF] Morgan | A homotopy for solving general polynomial systems that respects mhomogenous structures[END_REF] with an equation constraining the next point to rely at a given distance from the previous one. The prediction step is performed using a tangent predictor and the correction one uses a Newton-Raphson algorithm. More details about continuation techniques can be found in [START_REF] Allgower | Introduction to Numerical Continuation Methods[END_REF].

Bifurcation detection

In this paragraph, we will focus on turning points and branching points which are the main bifurcations encountered in this study. The detection of these bifurcation points is done by monitoring the determinant of the Jacobian matrix [J Y ] = ∂G ∂Y of (20): this determinant is null for such bifurcation points [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]. A Newton like algorithm, described in [START_REF] Nayfey | Applied nonlinear dynamics[END_REF], is used to determine accurately the solutions

(Y b , ω b ) such that G(Y b , ω b ) = 0 ∂G ∂Y (Y b , ω b )T = 0 T = 1 (21)
where T is the tangent direction.

The type of bifurcation is then determined by estimating the rank of the Jacobian matrix [J] = [J Y J ω ] (where [J ω ] = ∂G ∂ω ) which is of size ñ × (ñ + 1). This rank is exactly ñ for turning points and no further work has to be done since the arclength continuation scheme is able to handle such points. In the case of branching point bifurcations, this rank is at most ñ -1 and several linearly independent tangents T j respecting [J Y ]T j = 0 can be exhibited. Each of them corresponds to a prediction direction to explore. In this study, the prediction directions are computed using the eigenvectors Ψ i of [J Y ] associated with the zero eigenvalue at the bifurcation point [START_REF] Ribeiro | Non-linear free vibration of isotropic plates with internal resonance[END_REF].

Stability using Floquet theory

The stability of solutions is determined by using Floquet's theory. In this study, just the primary results of this method will be recalled; a more detailed presentation is available in [START_REF] Nayfey | Applied nonlinear dynamics[END_REF].

To do so, the first order formulation of the dynamic system ( 6) is considered:

U = f (U, t) (22)
with U a vector combining X and Ẋ.

If U * (t) is a periodic solution of ( 22) with period T * and z(t) is a small perturbation added to U * , the equation governing its evolution is derived from

(22): ż = ∂f ∂U (U * , t)z (23)
which is linear with periodic coefficients. Therefore, any solution z(t) corresponding to any initial perturbation is a linear combination of m (i) (t), solutions of (23) at time t having all its components but the i-th null at time t = 0: m (i)

j (0) = δ ij , 1 ≤ i, j ≤ 2n.
The growth of z is then directly related to the one of m (i) which is in turn evaluated via the eigenvalues of the monodromy matrix

[Mon] = m (1) (T * ), . . . , m (2n) (T * ) (24)
This was demonstrated and stated properly in the and on the type of bifurcation (through the way they cross the unit circle in complex plane) [START_REF] Nayfey | Applied nonlinear dynamics[END_REF][START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF].

Some other techniques [START_REF] Lazarus | A harmonicbased method for computing the stability of periodic solutions of dynamical systems[END_REF][START_REF] Laxalde | Étude d'amortisseurs non-linaires appliqus aux roues aubages et aux systmes multitages[END_REF] exploit the HBM results directly to determine the cycle stability. The authors confronted two methods. First a computation using Floquet basic theory relying on time integrations; the cycle U * is then rebuilt using the HBM coefficients. Second, the perturbation technique proposed by D. Laxalde in his PhD work [START_REF] Laxalde | Étude d'amortisseurs non-linaires appliqus aux roues aubages et aux systmes multitages[END_REF] was implemented. Each time, both methods gave the same results.

Application to backbone curves

The aim of this section is to build the NNMs of the undamped structure [START_REF] Morgan | A homotopy for solving general polynomial systems that respects mhomogenous structures[END_REF]. The results will be displayed mainly in the Frequency Energy Plot (FEP) which is one of the most appropriate ways to depict the frequency-energy dependence of such modes [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF]. First, a brief recall of the linear study results is proposed. Then, a global analysis step is performed and exhibits many solutions which are grouped with respect to their energy, creating energy branches. Each solution branch is finally followed and analyzed using tools of section 3.2. One of the branches and its bifurcations is analyzed in a detailed way before the drawing of the full bifurcation diagram.

Linear system

The linear system has the following pair of eigenvector and eigenvalues (eigenvectors normalized with respect to mass matrix), i.e. linear normal modes (LNM) :

ω 1 = 93.63 rad.s -1 , L 1 = {1, 1, 1, 1} t /2 ω 2 = 95.20 rad.s -1 , L 2 = {1, 0, -1, 0} t / √ 2 ω 2 = 95.20 rad.s -1 , L 3 = {0, -1, 0, 1} t / √ 2 ω 4 = 96.75 rad.s -1 , L 4 = {-1, 1, -1, 1} t /2 (25)
Using classical denominations [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF], one can state that • L 1 is a "mode 0" were all the coordinates have the same amplitude. This is a non-degenerate mode with no nodal diameter.

• L 2 and L 3 are degenerate "cos-1" and "sin-1" modes with 1 nodal diameter. According to [26, p. 202], the corresponding NNMs should have the same backbone (n = 4 × p × 1).

• L 4 is a "mode n/2" that exists because the system has an even number of blades. All coordinates vibrate with the same amplitudes but opposite phases, generating 2 nodal diameters.

These notations will be useful for backbone analysis in section 4.3.

Global analysis

A global analysis method as explained in 3.1 was applied to the free system [START_REF] Morgan | A homotopy for solving general polynomial systems that respects mhomogenous structures[END_REF] in order to build backbone curves. As the system is already in a polynomial form, no further work was done on equations. The total degree homotopy method was applied recursively for ω between 93.0 rad.s -1 and 100.0 rad.s -1 with a 0.1 rad.s -1 step. For each step 3 4 = 81 paths are followed. This computation took 758 s (about 12 min 38 s using an Intel Core 2 Duo E8400 (3Go RAM)) that is about 11 s per step and led to a great amount of solutions as depicted in Fig.

Besides the great number of solutions, one can

see that plots for each A i look similar which is in accordance with the cyclic symmetry of the structure.

Each solution returned by the global analysis algorithm was refined using a local Newton-Raphson method with a drastic 10 -12 criterion on residual;

no difference between raw GA results and refined ones appeared. We then ensured that each solution kept by the global analysis algorithm was different from the others from each step and counted them:

there are up to 81 solutions for pulsations higher than 99.4 rad.

-1 ! This justifies the use of the basic total degree homotopy method (same number of zeros for Q and P). This step also exhibited disconnected solutions. To have a better understanding of these numerous solutions, they were grouped with respect to their energy level. These energy levels are represented in a FEP, drawing branches that are named, as shown on Fig. 4. Each branch is named n.p 1 .p 2 . n denotes the number of nodal diameters of the main branch and a .p i extension is added each time a bifurcated branch occurs. d1 and d2 refers to the disconnected branches and NS to the null solution. Symbols on this figure refer to the number of solutions that have the same energy level for a given pulsation ω.

[Figure 3 about here.]

[Figure 4 about here.]

Branches analysis using continuation

Every solution corresponding to each energy branch has been continued and analyzed with regards to its stability and bifurcations using tools of section 3. Solutions underlying this energy branch 1'.1 will bifurcate once again for ω 8 = 99.36 rad.s -1 , that is for the same pulsation at which branches d1 and d2

arise.

[Figure 6 [Figure 8 about here.]

Conclusion

Table 1 and Fig. 9 summarize the results of the previous analysis carried out on all other branches. This simple structure exhibits in fact a very rich modal behavior: a great number of solutions branches can be captured using continuation and bifurcation analysis. Due to cyclic symmetry, they draw a less important but still great number of energy branches. In addition to these numerous branches, the use of a GA method lets us intercept two disconnected branches d1 and d2 that would have been missed otherwise.

[Table 1 about here.]

[Figure 9 about here.]

Application to a forced case

The previous section exhibited numerous NNMs which will lead to complicated structures of resonance in forced cases. The aim of this section is to illustrate this complexity through an example.

The system is therefore excited with a force taking the branch 1' shape, thus vector F e0 of Eq. ( 13) is defined by F e0 = a f {1, 1, -1, -1} t . Simulations are carried out for a f = 0.25 N.kg -1 .

First, the principal non-linear response has been computed by continuation. This solution, called Sol-A, corresponds to a motion with the 1' mode shape. The solution and its stability are represented in a FEP on Fig. 10.

By monitoring the determinant of the Jacobian matrix, 4 bifurcation points have been computed for Sol-A, they are also represented on Fig. The same representation code is used for the following figures (Fig. 11 to Fig. 13). One can then see on Fig. 10 that sol-A is unstable even for relatively low energy levels.

[Figure 10 B2 (or B3) leads to an unstable bifurcated response, Sol-C, positioned around the backbone curve of branch 1'.1.1 (Fig. 13), and B1 leads to an unstable bifurcated response, Sol-D, which seems to be positioned close to the backbone curve of branch 1' (Fig. 13).

[Figure 12 which seems reasonable with regard to the information it provides.

[Figure 13 about here.]

All the solutions are also depicted in a frequencyamplitude plot in Fig. 14 and Fig. 15. Two different figures were used for the sake of clarity. Only the amplitude of the first sector is considered.

[Figure 14 about here.]

[Figure 15 about here.]

In order to fully validate the stability analysis, direct numerical integrations have been performed.

Initial conditions were provided by HBM solutions and simulations were carried out over 1000 periods.

Results are plotted only for the five last periods.

Fig. 16 (resp. Fig. 17) shows the temporal evolution of both numerical integration and HBM solution for a stable point of Sol-B (resp. Sol-E) at ω = 99.53 (resp. ω = 98.84). A good agreement can be observed for both simulations, thus confirming the validity of the proposed approach.

[Figure 16 Applied to this case, global analysis has two virtues: first of all, finding the disconnected solutions, that is guarantying that all the possible stable states will be considered when dimensioning the structure.

The second one is not theoretical, but practical: if one can detect all the bifurcation points by carrying out a bifurcation analysis along the continuation process, it is easy to skip such points. For example, bifurcation points A2 and A4 being turning points with two additional tangent for each, the determinant of the Jacobian gets null but does not change sign. An automated script based on this criterion had therefore missed Sol-L2 and Sol-L3 branches and they were in fact found by the GA approach.

Conclusion

This paper applies both a classical continuation method and an original global analysis approach to a cyclic symmetric structure. These methods enable the drawing of a complete bifurcation diagram in a reasonable amount of time even in the case of numerous bifurcations and existence of disconnected solutions. Moreover it can intercept branches that would have been missed by the continuation and bifurcation analysis scheme. The forced case example also emphasizes the usefulness of NNMs to analyze a non-linear structure.

The ability of the proposed global analysis method to exhibit disconnected solutions being proved, further work now should focus on developments allowing its application to larger systems and to quasi-periodic solutions search.

A. Definitions of parameters in the equation of motion

This Appendix gives a symbolic and numeric expression for coefficients used in the equation of motion [START_REF] Huber | Polyhedral end games for polynomial continuation[END_REF].

m = ρ P Φ 2 dP = ρhL x L y 5 c = k m Φ(x r , y r ) 2 = 1 256 k m F ej = Φ(x f , y f ) f j (t) m = f j (t) m α = 5 3
Eh 2 ρL 4

x (1ν 2 )

β = 8E ρL 4 x (1 -ν 2 ) (26)
The damping coefficient δ is defined as:

δ = √ α 200 ( 27
)
Here are the numerical values used in this study:

the geometric and physical parameters are affected by the following set of values:

L x = 1.5 m, L y = 0.3 m, h = 0.03 m, x r = L x /4, y r = 0, x f = L x , y f = 0, E = 210 GPa, ν = 0.3, ρ = 7800 kg.m -3 , k = 8 • 10 5 N.m -1 (28)
which correspond to the following values for the different parameters in equation ( 6) : A 2 (mm)

α = 8.7662 • 10 3 s -2 , c = 148.36 s -2 , β = 4.6752 • 10 7 m -2 .s -2 (29)
A 3 (mm) A 4 (mm) ω (rad. -1 ) ω (rad. -1 ) ω (rad. -1 ) ω (rad. -1 )
Figure 3: GA applied to free system: solutions components

A 1 A 2 A 3 A 4 1' 1'.1 1'.1.1 Stable Unstable ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 2 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 3 ω 4 ω 4 ω 4 ω

  about here.] [Figure 2 about here.] The remainder of this article will be devoted to studying the case of a system composed of n b = 4 identical sectors, whose displacements are interpolated by a single Ritz function approximating the first bending mode thus yielding a non-linear problem with n = 4 dofs that will serve as an example for the present work. The Ritz function is Φ(x, y) = ( x Lx ) 2 (consistent with the clamping at

Floquet's

  theorem (1883). Eigenvalues of the monodromy matrix are usually complex and are called Floquet coefficients. They give us information on the stability of the solution (U * is unstable as soon as one of them has a magnitude greater than 1)

2 .

 2 As the previous global analysis step provided startpoints for each solution branch, the structure symmetry was not used to reduce the number of solutions to study; this would indeed have required a fine bifurcation analysis in conjunction with the symmetry consideration to determine the number of solutions appearing at each bifurcation point. As discussing each branch would be too long, one proposes to focus on branch 1' and its bifurcated branches. The same analysis was carried on for each branch and the results are summed up in subsection 4.4. For each energy branch 1', 1'.1 and 1'.1.1, one of the solutions matching this energy evolution is picked up as a representative of all the other solutions with the same energy evolution. Each time, a figure divided in four frames presents its HBM components evolution along with ω in frame (a); in frame (b) the evolution of its decomposition on LNM in percent is drawn: the solution components are projected on LNM of subsection 4.1 and these "modal coordinates" are brought back to percent by dividing by the sum of their absolute value and multiplying this by 100 . This will let us conclude on similar or non-similar property of the mode. Finally, Floquet coefficients are plotted in frames (c) and (d): frame (c) shows their norm along ω while frame (d) represents them in the complex plane superimposed on the unit circle. In order to get accurate Floquet coefficients, solutions initially computed using only the first harmonic A 1 were refined using A 1 , A 3 and A 5 components.Branch 1' matches 4 different solutions with the same energy that contain both L 2 and L 3 contributions with equal weight (Fig.5 (a)). These 4 solutions are identical with respect to a circular shift of tha A i components. Fig.5 (b)shows that it is a similar mode. Fig.5 (c) and (d)show that this mode bifurcates and become unstable for ω 3 = 95.98 rad.s -1 generating one new energy branch, namely branch 1'.1. In fact, each of the 4 solutions underlying energy branch 1' bifurcates and gives birth to two half branches that match the branch 1'.1 energy level.[Figure5about here.] Branch 1'.1 matches 8 solutions that come from branch 1' bifurcation. These solutions include a LNM L 1 contribution (Fig. 6 (b)) that becomes greater as ω increases. This indicates that this is not a similar mode. This branch is first unstable and becomes stable for ω 4 = 96.75 rad.s -1 , which is the frequency at which branch 2 arises. It then gives birth to a new branch (branch 1'.1.1) that "takes advantage" on the existence of this new energy level (of branch 2) as explained in the next paragraph. In the state space, each of the 8 half branches matching branch 1'.1 energy branch bifurcates once and gives birth to 2 new half branches.

  Fig. 8 recapitulates branch 1' bifurcations by showing the four initial solutions that match this energy branch and the subsequent bifurcated solution branches. It is obvious via this figure that each

Fig. 9

 9 Fig. 9 exhibits branches similar to Fig. 4 but this time, points are linked continuously thanks to the continuation step. Moreover, this figure provides information on stability and not on the number of solutions sharing the same energy branch anymore.

  10 labeled from A1 to A4. By studying the rank of the Jacobian matrix, one can tell that points A1 and A3 are branching points (rank([J]) = 7) and points A2 and A4 are not only turning points but also branching points (rank([J]) = 6). Stable parts of the energy branch are denoted by a thick (blue) line on Fig. 10 while unstable ones are represented by a thin (red) one. Underlying NNMs are plotted using dashed dot (black) lines. Squared points denote branching points while circled ones indicate Hopf bifurcations.

  about here.] Starting from A2 (or A4), the two bifurcated branches of solution have been computed (Sol-L2 and Sol-L3). The representation of these solutions in the amplitude frequency diagram is not very convenient because they superimpose themselves with a part of Sol-A, the only difference being that curves are not covered in the same way. To overcome this issue, we chose to represent these solutions in a FEP in Fig. 11. In this FEP, Sol-L2 and Sol-L3 are superimposed, but one can see that these solutions possess stable parts and are positioned around backbone curves of modes with 1 nodal diameter (branch 1), therefore one can conclude that these two solutions correspond to resonances around nonlinear modes with shape L2 and L3. [Figure 11 about here.] Starting from A1 (or A3), a bifurcated branch has been computed (Sol-B), and its representation in the FEP is given on Fig. 12. Sol-B is positioned around the backbone curve branch 1'.1. A stability analysis reveals that Sol-B possesses stable parts and also unstable parts generated through Hopf bifurcation. Several branching points have been detected for Sol-B, namely B1, B2 and B3 (Fig. 12).

  about here.] Finally, global analysis reveals another kind of solution, disconnected from the other. They are represented by closed curves in the FEP. This solution family, termed Sol-E, is positioned near the backbone branch 2.1 and is depicted on Fig. 13. Using formulation (11), 6561 paths are followed for each tested ω. This takes approximately 1 h 15 min per ω value using an Intel Core 2 Duo E8400 (3Go RAM)

  about here.] [Figure 17 about here.] This example exhibits a very complicated structure of resonance which is closely related to the backbone curves previously presented. The FEP appears to be a convenient tool for identifying nonlinear modes involved in non-linear resonances.

Figure 1 :Figure 2 : 6 A 1

 1261 Figure 1: Bladed wheel model used for establishing equations

Figure 4 : 2 ( 3 Figure 5 4 Figure 6 2 (

 4235462 Figure 4: GA for NNMs: energy branches

Figure 7 :

 7 Figure 7: Branch 1'.1.1 analysis

Figure 9 :- 1 )Figure 10 :- 1 )Figure 11 :- 1 )Figure 12 :- 1 )Figure 13 :Figure 14 :Figure 15 :Figure 16 :Figure 17 :

 911011111211314151617 Figure 9: Backbones FEP with stability

Table 1 :

 1 Backbones branches characterization sum up

	Branch Nsol ω i (rad.s -1 ) Sim/Non-sim
	0		ω 1 = 93.63	Similar
	1		ω 2 = 95.20	Similar
	1'		ω 2 = 95.20	Similar
	1'.1		ω 3 = 95.98	Non-similar
	1'.1.1	16	ω 4 = 96.75	Non similar
	2		ω 4 = 96.75	Similar
	2.1		ω 5 = 97.50	Non similar
	2.2		ω 5 = 97.50	Non similar
	2.3		ω 6 = 98.27	Non similar
	2.3.1		ω 7 = 98.45	Non similar
	d1		ω 8 = 99.36	Non similar
	d2		ω 8 = 99.36	Non similar

Preprint submitted to International Journal of Non-Linear MechanicsDecember 8, 2010

List of Figures