Samuel Quaegebeur
email: samuel.quaegebeur@ec-lyon.fr

Benjamin Chouvion

Fabrice Thouverez

Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties

Keywords: Cyclic symmetry, Friction nonlinearities, Traveling and standing wave excitations, Harmonic Balance Method

Displacement control of bladed-disks is of primary stake for turbo-engineers. Such structures present cyclic symmetric properties that allow some specific reduction techniques. For linear problems, the equation of motion projected on spectral components, also called nodal diameters, gives a system of equations in which the unknowns are uncoupled. Each nodal diameter can therefore be considered independently. However in real applications, the presence of nonlinear terms couples the different nodal diameters and makes the spectral equation of motion more complex to handle. This paper deals with this difficulty and presents two main results. It first gives an analytical derivation to determine which nodal diameters get coupled by friction nonlinearities. Such procedure reduces the size of the model but also the number of unknowns in the system by considering only the interacting nodal diameters. This method is general and allows to tackle a wide range of industrial problems. However this may not be sufficient for an efficient resolution of the nonlinear system since the nonlinear forces must first be evaluated in the physical domain, and, a priori, for all the sectors of the cyclic system before being computed in the spectral domain. The second main originality of the paper is the development of different strategies on this matter. One is analytical, valid for any excitation and is straightforward to implement, while the others are based on specific assumptions on the deformed shape but offer further reduction.

The new methodologies are validated on a simplified bladed-disk with different excitation forces and different friction's laws. They show very good accuracy and a substantial computation time reduction.

Introduction

This paper focuses on cyclic structures such as turbo-engines, with the assumption that they are perfectly tuned. For a linear system, the equation of motion of such structures are, most of the time, projected on their spectral components (also called nodal diameters) counterpart [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF][START_REF]NASTRAN cyclic symmetry capability[END_REF]. The spectral system is advantageous: as each nodal diameter is decoupled from each other, it greatly facilitates the solution procedure. In turbomachineries, one stage of bladed-disk is excited by the airflow going through the previous stages that are either rotating or fixed in the frame of the system studied. Therefore the excitation can either be a standing wave, a traveling wave, or a combination of both. This external force forms a specific wave shape due to the cyclic property of the structures. For the forced response of a linear system, only the nodal diameters equal to the wave numbers composing the force are excited and must be considered.

In real applications, nonlinear effects may arise in the system due to large displacement or friction between different bodies. These lead to two difficulties. Firstly, the nodal diameters get coupled by the nonlinear effects and the number of equations to solve is the same as in the initial physical system of equations, even after a spectral projection. This usually makes the system too large to be solved efficiently. Secondly, the nonlinear forces must be given in their spectral form in the solution procedure whereas they are generally known explicitly in the physical domain. Computing the spectral components of the nonlinearities usually requires to project the displacement in the physical domain, to compute the nonlinear forces and then to project them back into the spectral domain. This set of calculations is also time consuming.

Structural cyclicity combined with nonlinearities have interested the scientific community for a long time. The work of Vakakis [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF] presents many interesting insights on the behavior of such systems. The concept of nonlinear normal mode [START_REF] Rosenberg | The Normal Modes of Nonlinear n-Degree-of-Freedom Systems[END_REF] has been used in [START_REF] Georgiades | Modal Analysis of a Nonlinear Periodic Structure with Cyclic Symmetry[END_REF] to study the modal interaction of the cyclic structure. Using cyclic symmetric properties in the case of polynomial nonlinearities, analytical derivations were provided in [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF] to predict which nodal diameters interact with each other and to provide the expression of the nonlinear effects directly in the spectral domain. This method offers a large system reduction and enables to recover complex phenomena such as internal resonances. Under the specific case of a traveling wave excitation, Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] made the assumption that the displacements, and hence the nonlinear forces, have the same shape as the excitation. This assumption simplifies the system and still shows great accuracy for friction nonlinearities [START_REF] Petrov | Analysis of sensitivity and robustness of forced response for nonlinear dynamic structures[END_REF]. In the same paper, Petrov proposed an extension for the standing wave excitation case: the methodology was adapted by increasing the size of the sector of reference. For this "supersector", the standing wave excitation can be seen as a traveling wave excitation with the wave number equal to 0. This strategy is theoretically correct but is not practical for finite element models of industrial structures because one must still consider several sectors in the model and its associated size may be prohibitive. Moreover, in the case of a structure with polynomial nonlinearities, complex phenomena can be obtained (such as internal resonance [START_REF] Nayfeh | Nonlinear Oscillations[END_REF], multiple solutions [START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF] or localization [START_REF] Vakakis | A Multiple-Scales Analysis of Nonlinear, Localized Modes in a Cyclic Periodic System[END_REF]). Under these circumstances, the solution may not be following a traveling wave shape anymore and therefore the strategy of [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] can no longer be used.

The objective of the current work is to provide a new procedure enabling to reduce efficiently the size of the model and the number of unknowns in cyclic structures undergoing nonlinearities and possible complex excitation wave forms. The focus will be on friction nonlinearities only but some of the developments presented in the paper could be extended to other nonlinear effects. Analytical derivations are first proposed to determine which nodal diameters interact with each other. By this mean, the number of nodal diameters to consider in the spectral system of equations is reduced. However, as the nonlinear forces must generally still be evaluated in the physical domain, it may therefore not be efficient enough for a large number of sectors. An analytical procedure is then proposed to reduce the number of sectors for which the nonlinear forces must be evaluated. These analytical derivations give a new general method of reduction, free of any assumptions. Finally, to go even further in the reduction process and consider a smaller number of sectors in the nonlinear effects evaluation, this paper will also provide a new strategy based on an assumption on the displacement shape.

The paper is organized as follows, Section 2 recalls the properties of cyclic symmetry and provides the spectral system of equations for any nonlinearity. Analytical derivations are made in Section 3 to explicitly provide which nodal diameters are going to interact with each other due to the nonlinear terms. This result is supplemented with an analytical procedure to reduce the number of sectors for which the nonlinear forces must be evaluated. Section 4 proposes other strategies to reduce the size of the problem: the work of [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] will first be recalled before giving details on how it can be adapted with the previously proposed analytical development. Section 5 validates numerically the new methodologies for different friction law's and different excitation forces by comparison with simulations on a full sized model.

Spectral equation of motion for a nonlinear cyclic symmetric structure

The purpose of this section is to briefly derive the equations of motion of a cyclic symmetric system written in its spectral basis. These well known steps are used to properly introduce the notations. Complete derivations can be found in [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF][START_REF]NASTRAN cyclic symmetry capability[END_REF].

Let consider a cyclic structure composed of N sectors such as the one depicted in Figure 1 (where the sectors are represented separated from one another). The system is supposed perfectly tuned. Each sector 1 2 3 4 5 6 [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] 8 For the decoupled system (the cyclic boundaries are not assembled), the equations of motion of the whole structure are

M ü + C u + Ku + f c + f nl = f ext , (1)
where M = (I N ⊗ M 0 ), and similarly for C and K. ⊗ denotes the Kronecker product and I N the identity 65 matrix of size N. The vector u = u T 1 , ..., u T N T represents the displacements of the degrees of freedom (DOFs) of the entire structure, u j denotes the displacement for only the j-th sector ( j ∈ 1, N ). The term f c contains the forces that the sectors apply on each other. The term f nl represents nonlinear forces and f ext denotes external excitation forces. The cyclic boundaries are supposed to be free of nonlinear forces.

From the theory given in [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF], the physical displacement u can be transposed into its spectral compo-70 nents, such that

u = F ũ, (2)
where ũ = ũT 0 , ..., ũT N-1 is the full spectral displacement component vectors ( ũk represents the spectral displacement of order k) and F is the Fourier matrix given by

F = 1 √ N         1 1 • • • 1 1 e iα • • • e i(N-1)α . . . . . . . . . . . . 1 e i(N-1)α • • • e i(N-1) 2 α         ⊗ I N i +2N b , (3)
N i and N b are respectively the number of internal and cyclic boundary (right or left) DOFs. The scalar value α is the sector angle given by α = 2π N . As explained in [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF], only the quantities ( ũk

) k∈ 0,K , with K = N 2 if N is even and K = N-1
2 if N is odd, are required to retrieve the full displacement. The others components are their complex conjugates:

ũN-k = ūk for k ∈ 1, M , (4)
with

M = N 2 -1 if N is even and M = N-1 2 if N is odd. Substituting u by F ũ and projecting (1) on F gives FMF ü + FCF u + FKF ũ + F (f c + f nl ) = Ff ext . (5)
This expression can be simplified with the following property

N-1 ∑ p=0 e -ipnα e ipmα =      0 if n = m N if n = m, (6)
where m and n are integers (in our application these will represent nodal diameters numbers). It gives

M ü + C u + K ũ + F (f c + f nl ) = Ff ext . (7)
The displacement of each u j is decomposed into three parts: its internal DOFs (subscript i), its right interface DOFs (subscript r) and its left interface DOFs (subscript l). In order to assemble the structure, the following displacement continuity must be respected

u j,l = u j+1,r 1 ≤ j ≤ N, (8)
with u N+1 = u 1 .

These conditions can be expressed in the spectral domain, and one gets for each spectral component

ũk,l = e i kα ũk,r k ∈ 0, N -1 . (9)
For a finite element model (FEM), the DOFs are usually expressed in a local reference frame which supplements (9) with a rotation matrix. This will be omitted for simplification purpose and we consider that the DOFs displacements are all expressed in the same global frame. The information on cyclic boundaries [START_REF] Nayfeh | Nonlinear Oscillations[END_REF] are then used to reduce further the system [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF]. For each spectral order k, one has

     ũk,i ũk,r ũk,l      = Bk   ũk,i ũk,r   , where Bk =      I N i 0 0 I N b 0 e i kα I N b      (10)
By means of a block diagonal matrix B composed of Bk k∈ 0,N-1 , the new vector of solution becomes B ũ.

It will be noted ũ for simplicity in the rest of the paper. Equation ( 7) is premultiplied by BT to give:

M ü + C u + K ũ + BT F (f c + f nl ) = BT Ff ext , ( 11
)
where M is a block diagonal matrix which contains the matrices Mk k∈ 0,K with Mk = BT k M 0 Bk (and similarly for C and K). With Newton's third principle, the forces between sectors compensate each other out and assembling the different sectors gives BT Ff c = 0. The new system to be solved is given by:

M ü + C u + K ũ + BT Ff nl = BT Ff ext . (12)
In turbine engines, the excitation forces arise from the previous stages of bladed-disks and thus lead to both traveling wave (the force rotates along the structure) and standing wave (the force is fixed in the structure) excitations. These excitations are usually characterized with a given wave number. This wave number will be noted h ex (0 ≤ h ex ≤ K). Numerical examples will be provided in Section 5 for standing wave excitations and also for a combination of both traveling and standing waves. The cyclic boundaries are supposed to be free of external and nonlinear forces. Therefore BT does not impact the forcing term and the nonlinearity.

In a general manner, a mono-harmonic external force applied on the first sector can be written as

f ext,1 (t) = 1 2 f a e iωt + fa e -iωt , (13)
where f a is a complex vector and ω is the excitation frequency. For a sector j, the excitation force f ext, j can be computed directly with the force applied on the first sector f ext,1 (t). For a traveling wave excitation, their relationship is

f ext, j (t) = f ext,1 t - α ( j -1) h ex ω (14)
A standing wave can be seen as the sum of two waves traveling in opposite direction. Its expression can be given by

f ext, j (t) = 1 2 f ext,1 t - α ( j -1) h ex ω + f ext,1 t + α ( j -1) h ex ω = cos (( j -1) αh ex ) f ext,1 (t) (15)
Next, we study how a specific excitation impacts its spectral components. Let consider a spectral component q ∈ 0, K . For this component the term Ff ext is noted fext,q (similar notation is used for the nonlinear forces later on) and it can be calculated using [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF]. It gives:

fext,q = 1 √ N N ∑ j=1 e -iq ( j-1)α f ext, j . (16)
Employing the property of (6) for the traveling wave excitation yields,

fext,q = 1 2 √ N N ∑ j=1
e -iq ( j-1)α f a e iωt e -iα( j-1)h ex + fa e -iωt e iα( j-1)

h ex = 1 2 √ N N ∑ j=1
f a e iωt e -iα( j-1)(q+h ex ) + fa e -iωt e iα( j-1)(-

q+h ex ) =      0 if q = h ex √ N 2 fa e -iωt if q = h ex , . (17)
The derivations for the standing wave case are similar and give

100 fext,q =      0 if q = h ex √ N 4 f a e iωt + fa e -iωt if q = h ex , (18)
In both cases, the q-th spectral component of the external force, fext,q , is not zero if and only if q corresponds to the wave number of the external force. Notice that f nl in [START_REF] Vakakis | A Multiple-Scales Analysis of Nonlinear, Localized Modes in a Cyclic Periodic System[END_REF] depends on u (=F ũ) and/or its time derivative. Without this term, a linear problem is recovered and each spectral order k is decoupled in Equation [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF]. However in the presence of the nonlinear forces, all the nodal diameters may be coupled and thus the size of the problem is not reduced.

Reduction procedure for friction nonlinearities

The purpose of this section is to determine which nodal diameters are getting coupled by friction nonlinearities. To determine this coupling, the following work is based on an analytical development presented in [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF] that only handled polynomial nonlinearities. The dynamics involved with polynomial or friction nonlinearities is different as the former represents smooth functions while the latter is irregular. After introducing the main result of [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF], we will propose its extension to deal with friction nonlinearities.

3.1. Coupling of nodal diameters through the nonlinear forces 3.1.1. Brief recall of the cyclic nonlinear forces for polynomial nonlinearities Full details of this short derivation are provided in [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF] in which polynomial nonlinear forces were studied.

Those can be written as:

f nl, j = k nl u p j , (19)
with p an integer and k nl a scalar value. The displacement u j of sector j in ( 19) is substituted with its spectral component counterpart [START_REF]NASTRAN cyclic symmetry capability[END_REF]. Then the nonlinear forces are projected on their spectral components and the q-th spectral forces becomes:

fnl,q = 1 √ N N ∑ j=1 e -iq( j-1)α k nl √ N p p ∏ m=1 N-1 ∑ k=0 ũk e ik( j-1)α = 1 √ N p+1 k nl N-1 ∑ k 1 =0 N-1 ∑ k 2 =0
...

N-1 ∑ k p =0 ũk 1 ũk 2 ... ũk p N ∑ j=1 e i( j-1)α((k 1 +k 2 +•••+k p) -q) . ( 20
)
Based on relation [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF], fnl,q is different from 0 if and only if the set (k 1 , k 2 , ..., k p ) satisfies the following

relation p ∑ m=1 k m ≡ q mod N. ( 21
) Since 0 ≤ k m ≤ N -1, the relation (21) is equivalent to p ∑ m=1 k m = q + jN for j ∈ 0, p -1 . ( 22
)
If this condition is satisfied, the spectral components (k m ) m∈ 1,p are coupled by the nonlinear forces and excite the nodal diameter q.

Application to friction nonlinearities

In mechanical engineering applications, friction forces are the nonlinear forces that two structural bodies in contact exert on each other. One common macroscopic way to describe the friction effects is to use

Coulomb's law      f nl,T < µ f nl,N if ur,T = 0 f nl,T = -µ f nl,N ur,T ur,T if ur,T > 0, (23)
where f nl,T is the tangential friction nonlinear forces, f nl,N is the normal forces perpendicular to the contact area, µ is the friction coefficient and ur,T is the relative tangential velocity between the two bodies in contact.

As illustrated in Figure 2, this law is highly irregular. Different strategies exist to smoothen the nonlinearity. For instance, Firrone and Zucca [START_REF] Firrone | Modelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks[END_REF] presented different models for one and two dimensional 125

contacts. An elastic contact was used in [START_REF] Firrone | The effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method[END_REF][START_REF] Mitra | Adaptive Microslip Projection for Reduction of Frictional and Contact Nonlinearities in Shrouded Blisks[END_REF] and leads to stick slip motion. The following regularized

Coulomb's law also proves to be interesting [START_REF] Laxalde | Non-Linear Modal Analysis for Bladed Disks With Friction Contact Interfaces[END_REF][START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF] for numerical simulations,

f nl,T = µ f nl,N tanh ur,T ε , (24)
where ε is a positive parameter which controls the slope near the zero relative velocity (see Figure 2). The smaller ε, the more the regularized function fits Coulomb's law. The hyperbolic tangent friction law is adopted in the rest of this section for analytical purpose. The 130 Taylor expansion of the hyperbolic tangent function is

tanh (z) = ∞ ∑ n=1 q n z 2n-1 , (25)
where q n is the Taylor expansion coefficients of order n. Using this Taylor expansion within [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF], and, in a similar manner as performed in [START_REF] Lemoine | Influence of fretting wear on bladed disks dynamic analysis[END_REF], projecting the nonlinear forces on the spectral component q (q ∈ 0, K )

gives:

fnl,q = µ f nl,N √ N ∞ ∑ n=1   q n √ Nε 2n-1 N-1 ∑ k 1 =0   N-1 ∑ k 3 =0   ... N-1 ∑ k 2n-1 =0   uk 1 uk 3 ... uk 2n-1 N ∑ j=1 e i( j-1)α((k 1 +k 3 +•••+k 2n-1 )-q)         . (26)
Therefore if there exists an n ∈ N * such that 2n-1

∑ m=1 k m ≡ q mod N, (27)
then fnl,q is different from 0 and directly excites the spectral component ũq . This Equation ( 27) provides all the coupling terms that the nonlinear forces can create and helps to determine which spectral components will be excited.

These analytical derivations have been obtained with a regularized Coulomb's law. However taking ε close to 0 results in the asymptotic case and remains valid with our approach. As it is important to generalize this result, numerical simulations that employ the Dynamic Lagrangian Frequency Time (DLFT) algorithm [START_REF] Nacivet | A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF] which assumes no regularization of the friction law will be performed in Section 5.

Overall, the methodology allows to reduce the system (12) by considering only the interacting nodal diameters. In practice, the Fourier matrix defined in (3) is reduced in the rectangular matrix F r that contains only the coupled diameters. The displacement of the structure is then found with

u = F r ũr + c.c, ( 28
)
where ũr is the spectral components of the displacements for the interacting nodal diameter up to K. c.c is the complex conjugate terms to account for (4). Therefore the system (12) becomes:

Mr ür + Cr ur + Kr ũr + FT r f nl = FT r f ext . (29)
Notice that ( 12) was obtained with the projection on the complete Fourier matrix (square and symmetric), whereas system (29) is obtained with the projection on the reduced Fourier matrix FT r which is no longer square. This single reduction may not be efficient enough due to the process required to evaluate the nonlinear forces. In numerical applications, one would first need to compute the physical displacement for all sectors, then apply the nonlinear law and finally project the force on the interesting nodal diameters. This may become prohibitive for a finite-element model of a real bladed-disk. Therefore Section 3.2 will address this issue and present a general reduction valid for any wave excitation.

Examples and algorithm

Let consider an external force with a wave number h ex . Without nonlinear coupling, only the spectral component ũh ex responds (see Section 2). k m can then only take the value h ex , and condition ( 27) is simplified in the following: if there is a n ∈ N * such that

(2n -1) h ex ≡ q mod N, (30)
then fnl,q is different from 0 and the spectral component ũq is excited. Although n is an integer (and can thus take an infinite number of values), the remainder of the Euclidean division of (2n -1) h ex by N is finite and all of its values are obtained for a finite number of possible n-values. When n is larger than a particular integer, called M 1 , then the same remainders are obtained once again. The value of M 1 can be obtained with the equation:

(2 (n + M 1 ) -1) h ex ≡ q mod N, (31)
which gives

2M 1 h ex ≡ 0 mod N. ( 32
)
The scalar value M 1 can then be evaluated easily. By varying n ∈ 1, M 1 we obtain all the possible values of the remainder of (2n -1) h ex by N for n ∈ N * . These values correspond to the interacting nodal diameters.

A numerical example will be provided later in this section.

Through this first iteration, we have obtained a set of nodal diameters {q 1 , • • • q P } that are coupled through an excitation h ex . One could then think that the following iterative procedure must be employed:

consider the new set of nodal diameters k m ∈ {q 1 , • • • , q P }, evaluate all the possibilities from ( 27), find new interacting nodal diameters, and so on, until convergence is reached and no more new nodal diameters are found with the algorithm. However such an iterative process is not needed: the nodal diameters found with the first step are necessarily equal to an odd multiplicative factor of h ex (with a possible shift of N): [START_REF] Poudou | Hybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems with Dry Friction Damping[END_REF] with the new possible values of k m . We split this general sum into multiple sums for each different nodal diameters. It gives

q i ≡ (2n i -1) h ex mod N, i ∈ 1, P . Consider Equation
m 1 ∑ q 1 + • • • + m P ∑ q P ≡ q mod N, (33)
with ∑ i m i = 2n -1. As explained previously each nodal diameter q is in fact an odd multiplicative factor times h ex . Thus (33) can be written as

m 1 ∑ (2n 1 -1) h ex + • • • + m P ∑ (2n P -1) h ex ≡ q mod N, (34)
Factoring out the value h ex yields

h ex P ∑ i=1 (m i (2n i -1)) ≡ q mod N, h ex 2 P ∑ i=1 m i n i - P ∑ i=1 (m i ) ≡ q mod N, (35)
since ∑ i m i = 2n -1, the last equation is just the remainder of the Euclidean division by N of an odd multiplicative factor times h ex . This equation is therefore equivalent to [START_REF] Poudou | Hybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems with Dry Friction Damping[END_REF], and thus solving it once enables to recover all interacting nodal diameters. The methodology can be generalized to multiple excitations with multiple wave numbers (see Appendix A).

Numerical example

Let consider the case N = 24 and h ex = 3. The computation of (32) gives M 1 = 4. The values of the term (2n -1) h ex for n ∈ 1, 4 gives [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Mitra | Adaptive Microslip Projection for Reduction of Frictional and Contact Nonlinearities in Shrouded Blisks[END_REF][START_REF] Gross | Effect of creep on the nonlinear vibration characteristics of blades with interlocked shrouds[END_REF]. As ũ15 = ū9 and ũ21 = ū3 (see Equation ( 4)), only the 3 rd and 9 th nodal diameters will get coupled. Then k m ∈ {3, 9}. Let detail (33) to justify that there is no need to re-apply the algorithm with this new set of values for k m :

m 1 ∑ 3 + m 2 ∑ 9 ≡ q mod N, m 1 ∑ 3 + m 2 ∑ (3 × 3) ≡ q mod N, 3 (m 1 + 3m 2 ) ≡ q mod N, 3 (2 (m 1 + 2m 2 ) -(m 1 + m 2 )) ≡ q mod N, (36)
The term m 1 + m 2 is equal to 2n -1 and the term 2

(m 1 + 2m 2 ) is even. As a consequence 2 (m 1 + 2m 2 ) - (m 1 + m 2 ) is an odd number. Computing the remainders of the division of 3 (2 (m 1 + 2m 2 ) -(m 1 + m 2 )) by
N will provide the same values as the first iteration.

Reduction of the nonlinear evaluation for a general wave excitation (Method 1)

The purpose of this section is to propose a general and efficient method to compute the term FT r f nl of system (29). This term represents the projection on a reduced Fourier matrix of the physical nonlinear forces for all sectors. Such operation is time consuming and is to be avoided as much as possible. The following proposes an analytical approach to reduce the number of sector for which the nonlinear forces must be evaluated.

Similarly to the displacement reduction in [START_REF] Charleux | Etude des effets de la friction en pied d'aube sur la dynamique des roues aubagées[END_REF], we have for the nonlinear forces,

f nl = F r fnl,r + c.c, (37)
where fnl,r is the spectral component of the nonlinear forces for the interacting nodal diameters. This equation can be split into

f nl =   f nl,m f nl,s   = 2   Re (F r,m ) -Im (F r,m ) Re (F r,s ) -Im (F r,s )     Re fnl,r Im fnl,r   , (38)
where f nl,m are the nonlinear forces applied on the sectors 1, N diam , with N diam the number of interacting nodal diameters pondered by a coefficient 2 if the nodal diameter is degenerated. f nl,s are the nonlinear forces applied on the remaining sectors N diam + 1, N . The spectral nonlinear forces are split into their real and imaginary parts denoted with the operators Re and Im. Notice that for non-degenerated nodal diameters, the imaginary part of the spectral nonlinear forces is equal to 0 and is thus removed from the equation. The term

[Re (F r,m ) -Im (F r,m )
] is a square matrix of size N diam × N diam . Inversing this matrix and substituting the first line of (38) into the second line gives:

FT r f nl = FT r   I N diam [Re (F r,s ) -Im (F r,s )] [Re (F r,m ) -Im (F r,m )] -1   P f nl,m (39)
Finally, following this substitution procedure, the system (29) is reduced in

Mr ür + Cr ur + Kr ũr + FT r Pf nl,m = FT r f ext . (40)
In contrast of system (29), this new system (40) does not require an evaluation of the nonlinear forces for the entire structure but only for the subset of sectors 1, N diam . Notice that even though the evaluation must be performed for several sectors, only the spectral matrices of the reference sector are needed.

The previous development can be extended to the entire set of interacting nodal diameters of h ex . The interacting nodal diameters are equal to a product of an odd number by h ex whith a possible shift of N (as shown in Section 3.1). It corresponds to the same wave moving faster and thus the reduced number of sectors N diam remains the same. Therefore the new method can also handle excitation forces which combine multiple traveling and standing wave excitations whose wave numbers are equal to an odd number times h ex modulo N.

This reduction is general, valid for any wave excitation and straightforward to implement in a numerical code. The choice was made in this paper to study nonlinearities coming from frictional effects. However, the results remain valid for any polynomial nonlinearity.

Other strategies will be presented in Section 4. But those require specific approximations to reduce further the computation of the nonlinear forces. Section 4 will first explain the strategy developped by Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] for both a traveling wave and a standing wave excitations. Then we will show how this strategy can be extended with the previous development that predicts the interacting nodal diameters. Moreover, considering an additional particular assumption on the deformed shape, an original method will be given to decrease further the size of the problem for the case of a standing wave excitation.

Numerical example

Consider N = 24 sectors and an excitation force which combines traveling and standing wave excitations with h ex ∈ {3, 9}. Only the third and ninth nodal diameters interact. Both of them are degenerated nodal diameters and thus N diam = 4. The matrix P of (39) is

P =      1 1 1      ⊗                     1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 -1 0 0 0 0 -1 0 0 0 0 -1 0 0 0 0 -1                     . ( 41
)
For better readability, the matrix is represented as a Kronecker product. The size of the first vector 3 × 1 comes from the h ex = 3 excitation.

Strategies for specific wave excitation

Before proposing two original methods of reduction, this section will first recall the strategies proposed by Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] for the case of standing and traveling wave excitations.

4.1. Summary of the article by Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] 4.1.1. Reduction for a traveling wave excitation

In the specific case of a traveling wave excitation, Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] developed a methodology to reduce the size of the problem [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF]. This methodology has been widely used, for instance to study shrouds [START_REF] Siewert | Multiharmonic Forced Response Analysis of a Turbine Blading Coupled by Nonlinear Contact Forces[END_REF], fretting of bladed-disks [START_REF] Lemoine | Influence of fretting wear on bladed disks dynamic analysis[END_REF] or creep effect [START_REF] Gross | Effect of creep on the nonlinear vibration characteristics of blades with interlocked shrouds[END_REF]. As opposed to the analytical development proposed in Section 3.1, the work in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] relies on the specific assumption that the displacement (and thus the nonlinear forces) has the same shape as the traveling wave excitation [START_REF] Firrone | The effect of underplatform dampers on the forced response of bladed disks by a coupled static/dynamic harmonic balance method[END_REF]. It uses the following assumption:

u j (t) = u 1 t - α ( j -1) h ex ω , (42)
which is fundamentally different from the proposed approach in Section 3. Both approaches are actually complementary, as it will be shown at the end of this section. To better understand the differences, the method of [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] is explained first.

To take into account the nonlinear effects, the Harmonic Balance Method (HBM) [START_REF] Kerschen | Nonlinear normal modes, Part I: A useful framework for the structural dynamicist[END_REF] is used, and, the steady-state solution on the first sector, supposed periodic, is sought as a Fourier series

u 1 (t) = c 0 + 1 2 N h ∑ n=1 c n e inωt + cn e -inωt , (43)
where N h is the total number of harmonics retained and c n are the Fourier components (notice that c 0 is real).

Due to the particular shape of the solution (42), the spectral displacement can be easily computed. Using the same derivations as the traveling excitation force [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF], we have for each harmonic n ∈ 0, N h and a given wave number h ex :

ũq = 1 2 √ N N ∑ j=1
c n e inωt e -iα( j-1)(q+nh ex ) + cn e -inωt e iα( j-1)(-

q+nh ex ) =                            √ N
2 c n e inωt + cn e -inωt if q = 0 and nh ex ≡ 0 mod N,

√ N 2 c n e inωt + cn e -inωt if q = N 2 and (nh ex ) ≡ N 2 mod N √ N
2 cn e -inωt if (q -nh ex ) ≡ 0 mod N, and q / ∈ 0,

N 2 √ N 2 c n e inωt if (q + nh ex ) ≡ 0 mod N, and q / ∈ 0, N 2 0 otherwise . ( 44
)
In practice the user provides the number of harmonics considered (N h ) and evaluates for each harmonic 240 the remainder of (nh ex ) by N to determine q. For each harmonic n, a single value of q is obtained. This process must be applied for all n ∈ 0, N h . In the end, a table of correspondence can be created between the harmonics number and the nodal diameter. For each harmonic n, we finally have

ũq =      √ N 2 cn e -inωt or √ N 2 c n e inωt if q ∈ 1, M √ N 2 c n e inωt + cn e -inωt if q ∈ 0, N 2 . ( 45
)
This equation is substituted in [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF] and the system is projected on the exponential basis e inωt n∈ 0,N h with the scalar product:

f , g = 1 T T 0 f (t) ḡ (t) dt, ( 46
)
where T is the fundamental time period equal to 2π ω . Finally, the system (12) which initially required the computation of all nodal diameters, each of whom was expanded up to N h harmonics, is now simplified to 245 (inω) 2 Mq -(inω) Cq + Kq √ N 2 cn + fnl,q , e -inωt = fext,q , e -inωt ∀n ∈ 0, N h .

This means that each harmonic only affects one nodal diameter. The harmonic coefficients c f nl ,n of the nonlinear force on the first sector f nl,1 are evaluated through an AFT procedure [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF]. It gives

f nl,1 (t) = c f nl ,0 + 1 2 N h ∑ n=1 c f nl ,n e inωt + c f nl ,n e -inωt . (48)
The nonlinear forces for the remaining sectors are obtained with the assumption (42)

f nl, j (t) = f nl ( u j (t)) = f nl u1 t - α ( j -1) h ex ω = f nl,1 t - α ( j -1) h ex ω . (49)
Finally fnl,q , e -inωt is obtained with a similar development as the one used for ũ (Equation ( 43) to (45)),

fnl,q =      √ N 2 c f nl ,n e -inωt or √ N 2 c f nl ,n e inωt if q ∈ 1, M √ N 2 c f nl ,n e inωt + c f nl ,n e -inωt if q ∈ 0, N 2 . ( 50
)
This reduction proposed by Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] is powerful: instead of having a full system to solve of size N dof × N × (1 + 2N h ) (with N dof the number of DOFs in a single sector), we have a reduced system of size N dof ×

(1 + 2N h ). It offers a reduction by N. This reduction could have been expected since the assumption of a traveling wave displacement implies that knowing the displacement on one sector provides the displacement for the full structure. However, seeking the solution as a traveling wave shape is detrimental for bifurcation analysis where symmetry-breaking bifurcation points destroys the initial shape of the solution [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. Moreover in the case of a standing wave excitation the inital hypothesis is broken and the traveling assumption is, a priori, not valid.

Numerical example

As a numerical example, consider the case of N = 24, h ex = 3 and N h = 7 harmonics. The reduction process of matching one harmonic to a nodal diameter is depicted in Table 1. The harmonics n ∈ {1, 7} activate the third nodal diameter (see Equation ( 44)), harmonics {3, 5} activate the ninth nodal diameter and so on.

Harmonic number n 1 2 3 4 5 6 7

Nodal diameter q 3 6 9 12 9 6 3

Table 1: Table of correspondence between the harmonic number and the nodal diameter with a traveling wave assumption.

Reduction for a standing wave excitation

In order to handle standing wave excitations, Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] proposed an alternative to apply its strategy.

The N initial sectors are decomposed into M groups of N 1 sectors (N 1 depends on h ex and N). The system remains cyclic. The standing wave excitation with an initial wave number h ex applied on the N sectors is transformed into a traveling wave with zero nodal excitation applied on the M sectors. The previous strategy explained in Section 4.1.1 can then be employed for M sectors with h ex = 0. This method is interesting but 265 difficult to use in practice as one must create a supersector composed of N 1 sectors and the final size of the system remains large in the case of a finite element model of an industrial structure. In order to illustrate this methodology, Figure 3 represents a standing wave excitation with h ex = 3 and N = 24. We can see that if we consider M = 3 groups of N 1 = 8 sectors then this "supersector" undergoes a traveling wave with h ex = 0.

In order to compute the spectral nonlinear forces, the nonlinear forces must be computed for N 1 sectors in

Prediction of nodal diameters interaction combined with Petrov's method [7] (Method 2)

The purpose of this section is to combine the work of Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] presented in Section 4.1.1 with the possibility to predict which nodal diameters will interact (see Section 3.1).

The work in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] detailed in 4.1.1 for a traveling wave excitation enables to consider only specific harmonics in the HBM for a set of nodal diameters. The method in Section 3.1 enables to reduce the number of nodal diameters to consider. Combining them both lowers the overall number of unknowns and allows to evaluate the nonlinear forces on a single sector. Based on the previous example (N = 24 and h ex = 3), the algorithm of Section 3.1 provided the interacting nodal diameters {3, 9} (see Section 3.1.3). In Table 1, nodal diameters 6 and 12 can thus be removed from the reduction procedure. Notice that only odd number of harmonics are thus retained. This result can be explained by the fact that since the friction is antisymmetric only odd harmonics will respond on the main branch of solution of the system. Although the reduction is very interesting, it is only true for a traveling wave excitation and assumes that the dynamic of the structure shows a traveling shape under forcing.

In the case of a standing wave excitation, Section 3.1 allows to determine exactly which nodal diameters will respond and thus there is no need to consider a "supersector" (the description of one single sector is sufficient). This enables to reduce drastically the size of the model. Based on the previous example (N = 24 and h ex = 3), the method in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] required a supersector composed of N 1 = 8 standard sectors. The unknowns of the problem are then the displacements of these N 1 sectors. In the new method, only the reference sector must be considered and the number of unknowns is reduced by 2 (only ũ3 and ũ9 must be evaluated).

However the evaluation of the nonlinear force must still be performed on N 1 = 8 sectors and this may be time consuming. The next section proposes a new strategy to deal with this issue.

Reduction of the nonlinear force evaluation in the case of a standing wave excitation (Method 3)

The purpose of this section is to propose a new reduction of the number of sectors for which the nonlinear forces must be evaluated to solve [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF].

Consider a structure of N sectors excited with a standing wave force [START_REF] Mitra | Adaptive Microslip Projection for Reduction of Frictional and Contact Nonlinearities in Shrouded Blisks[END_REF]. The spatial distribution is embedded in the term cos (( j -1) αh ex ). An example of a standing wave excitation is given in Figure 3. The new proposed strategy is based on symmetric properties that allow to determine a pattern in the different levels of the standing wave excitation. By this mean, it only requires the computation of a small subset of sectors to retrieve the solution for the entire structure.

The first property that offers a standing wave shape is spatial periodicity whose period is composed of

N 1 sectors. This translates into cos (( j -1) αh ex ) = cos (( j -1 + N 1 ) αh ex ) ⇐⇒ N 1 αh ex ≡ 0 mod N (51)
The spatial periodicity N 1 corresponds to the minimum value for which Equation ( 51) is satisfied. In Figure 3 we have represented this condition with a red dashed square. This condition is equivalent to the one provided by Petrov (M group of N 1 sectors for a zero nodal diameter excitation). Within a group of N 1 sectors, the spatial distribution of the excitation is cos (2π ( j -1) /N 1 ). Depending on the value of N 1 , one considers three different cases: "N 1 is odd", "N 1 is even" and "N 1 is divisible by 4". For each of these cases, different 305 reductions based on symmetric properties are possible. When N 1 is odd, see Figure 4a,

N 2 = (N 1 + 1) /2
sectors are sufficient to describe the entire force distribution (through the use of an horizontal symmetry).

When N 1 is even and not divisible by 4, see Figure 4b only N 2 = N 1 2 + 1 /2 sectors are required to determine the entire spatial distribution. Finally, when N 1 is divisible by 4, see Figure 4c, only N 2 = N 1 /4 sectors need to be computed (through the use of an horizontal symmetry followed by a vertical one). For this last 310 case there is no need to consider the sector situated on the y-axis since cos (( j -1) αh ex ) = 0 for this sector and thus the excitation is zero.

In the particular case that N = 24 and h ex = 3, such as illustrated in Figure 3, N 1 = 8. It corresponds to Figure 4c and therefore N 2 = 2 as represented by the green dashed box in Figure 3. As a consequence, the standing wave excitation presents only N 2 non-zero different amplitudes. Shifting 315 these N 2 forces properly, one can retrieve the forces on the N 1 sectors. The pattern then repeats itself for the M groups (the structure was decomposed into M identical "supersectors" of N 1 sectors).

Let denote

f ext = f T ext,1 , • • • , f T ext, j , • • • , f T ext,N
T the standing wave excitation for all sectors. It can then be written as

f ext = Pf ext,r , (52)
where

f ext,r = f T ext,1 , • • • , f T ext,N 2
T is the reduced excitation forces and P is the pattern matrix (determined with the aforementioned symmetric properties) of size N × N 2 . In the numerical example given before, P is obtained with

P =      1 1 1      ⊗                     1 0 0 1 0 0 0 -1 -1 0 0 -1 0 0 0 1                     P m . ( 53
)
and is composed of the matrix P m that defines the symmetric properties for the N 1 = 8 first sectors, and a Kronecker product (of size M = 3) to repeat the pattern M times. For finite-element models, this matrix must be completed with rotation matrices that express local reference frames between the different sectors.

One now makes the assumption that the displacement and the nonlinear forces have a shape similar as the external force (52):

u = Pu r , (54)
with u r = u T 1 , • • • , u T N 2
T the reduced displacement. Combining this assumption with the result of Section 3.1 (see Equation ( 28)) yields ũr = Fr Pu r .

(55)

Finally the system (29) is once more reduced and gives (40) with the new matrix P. The nonlinear forces must only be evaluated for N 2 sectors. This methodology, with the assumption (55), will be called "Method 3" in the remaining of the paper. It simplifies the AFT procedure (explained in Section 4.1.1) by reducing the number of sectors that must be considered.

Performance of the different methodologies

The purpose of this section is to compare the reduction performance of the reference solution (modeling of the full system), the method in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF], and the proposed methodologies. The HBM is employed in all solution procedures. Section 5 will provide numerical comparison on frequency response functions. Tables 2 and3 summarize the advantages of these different methods in the cases of traveling and standing wave excitations.

A ranking based on numerical performance is also provided.

For the traveling wave excitation case (see Table 2), the method in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] with the HBM enables to associate each N h harmonic to a single nodal diameter. Method 2 uses the same results but only takes into account the harmonics which are associated to interacting nodal diameters (determined by Section 3.1). This reduces the number of harmonics that needs to be considered (size N h,2 < N h ). Both of these methods are more efficient than Method 1 due to the strong assumption that the displacement has a traveling wave shape.

For the standing wave excitation case (see Table 3), the new methodologies (Methods 1, 2 and 3) present many interesting advantages compared to [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF]. Firstly, they only require to model a single sector before applying the cyclic symmetry procedure (whereas the former needed to model N 1 sectors). This corresponds to a reduction by N 1 of the model. Secondly, the size of the problem to solve is also greatly reduced.

For the former method, a solution for the full supersector was required. With the new methods, its size is reduced to N diam (smaller than N 1 ) nodal diameters. The number of sectors for which the nonlinearities must be evaluated depend on the chosen method. With the assumption of a particular shape of the solution, Method 3 requires an evaluation of the nonlinear forces for N 2 (smaller than N diam ) sectors. Method 1, is the second most efficient method and requires the evaluation for N diam sectors. Lastly, comes Method 2 with an evaluation for N 1 sectors.

Notice, that in the case of an excitation with traveling and standing waves combined, only Method 1 is correct and it keeps the same reduction capabilities as the ones given in Table 2 and3.

Methodology

Size of the reference sector Nonlinear solver (HBM) Rank Number of unknowns AFT procedure

Reference N × N dof (2N h + 1) × N × N dof Evaluation for N sector 4 Former Method [7] N dof (2N h + 1) × N dof Evaluation for 1 sector 2 Method 1 N dof (2N h + 1) × N diam × N dof Evaluation for N diam sectors 3 Method 2 N dof 2N h,2 + 1 × N dof Evaluation for 1 sector 1
N × N dof (2N h + 1) × N × N dof Evaluation for N sector 5 Former Method [7] N 1 × N dof (2N h + 1) × N 1 × N dof Evaluation for N 1 sectors 4 Method 1 N dof (2N h + 1) × N diam × N dof Evaluation for N diam sectors 2 Method 2 N dof (2N h + 1) × N diam × N dof Evaluation for N 1 sectors 3 Method 3 N dof (2N h + 1) × N diam × N dof Evaluation for N 2 sectors 1
Table 3: Comparison of the different methods for a standing wave excitation.

To better assess the performances of the new methods for a standing wave, the numerical values of N 1 , N 2 and N diam are provided in Table 4 for different values of N and h ex . For N = 24 and h ex = 2 the former method [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] needed to create a supersector of 12 sectors and then to compute the nonlinear forces for these The Method 1 (nonlinear evaluation on the N diam sectors) is general and can handle a wide range of industrial problems while reducing significantly the size of the model. This methodology is expected to be able to recover internal resonance and follows bifurcated branches from symmetry-breaking bifurcation 360 points [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. The reduction proposed in [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF], valid in the specific case of polynomial nonlinearities, is expected to be more efficient as the nonlinear forces are directly evaluated in the spectral domain. In Method 2

(respectively Method 3), the number of sectors for which the nonlinear forces must be evaluated is greatly reduced for a traveling wave excitation (respectively standing wave excitation). However the hypothesis made on the deformed shape of the structure may not be valid for complex behaviour (combined excitation, gyroscopic terms, ...).

Numerical examples

The previous Sections 3 and 4 have established three new methods for the reduction of models undergoing friction nonlinearities. Method 1 is general and makes no assumption. Method 2 combines predicting spectral component interaction with Petrov [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] methodology. Method 3 proposes an additional reduction of the number of sectors for which the nonlinear forces must be evaluated but requires the assumption that both the displacement and the external force have a standing wave shape. The purpose of this section is to compare the accuracy and reduction performances of these methods with respect to a reference solution (full solution of the entire system).

The model considered in this study is a simplified bladed-disk represented in Figure 5. This model has been previously used in [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF] to assess the performance of another reduction method based on a nonlinear superelement: the Component Nonlinear Complex Mode Synthesis (CNCMS) procedure. For the numerical values of the masses, spring stiffnesses and dampers, the reader is referred to Table 2 of [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF]. The mass values were chosen to be representative of a bladed-disk mass distribution. Spring stiffnesses values were calculated such that the first bending mode family matches the one of a linear bladed-disk, and modal damping of 0.1% was used, which is representative of classical structural damping. In the simulation the term f nl,N is taken equal to 1000 N and µ to 0.

Frictional effects modeled with a regularized Coulomb's law

Methods 1, 2 and 3 are first employed to simulate the structure illustrated in Figure 5, under a standing wave excitation of wave number h ex = 3. Its amplitude varies between 1 N and 5 N with a 1 N step. The frictional effects are modeled with a regularized Coulomb's law (defined in [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF]). For all methods, 15 harmonics are retained in the HBM procedure. This number is sufficient to reach convergence on the nonlinear effects evaluation. The results of the different simulations are illustrated in Figure 6. For both Petrov's method and the new methods, the responses match perfectly with the reference solution. The wide range of excitation amplitude allows to cover both small and high levels of nonlinearities. For the ũ3 component represented in Figure 6, the different peaks vary between 207 Hz (corresponding to a configuration obtained with low excitation values) and 205 Hz (slip case at high excitations). Even for a small excitation, a coupling occurs since the ninth nodal diameter responds even though only the third nodal diameter is directly excited.

Table 5 provides the number of unknowns and the computation time for the different methods. For the proposed methodologies (Methods 1, 2 and 3), the size of the model is reduced by 87.5% compared to [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF], which offers tremendous computing gain and is especially interesting for industrial FEM. Moreover the size of the problem to solve is divided by 2, which allows Method 1 to decrease the computation time by about 53%. For Method 3 in which the nonlinear forces are evaluated only for two sectors, this gain increases and reaches 65%. The only difference between Methods 1, 2 and 3 is the evaluation of the nonlinear forces and thus one can clearly see the negative impact of computing the nonlinear forces on too many sectors. The simulated dynamical response of the structure is compared in Figure 7 with another reduction procedure: the CNCMS procedure [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF]. The main idea of the CNCMS is to create nonlinear Craig-Bampton superelements (with the computation of nonlinear normal modes with fixed-boundaries). This reduction methodology is efficient as the nonlinear forces are evaluated only once for the nonlinear modes and are then substituted in the synthesis procedure [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF]. It enables to reduce drastically the computation time as shown in Table 5. However these nonlinear modes with fixed boundaries may not correctly capture the dynamics of the system when the nonlinearities are important or when the vibratory energy of the disk is large (these limitations were underlined in Section 5.3 of [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF]). The CNCMS manages to capture coupling between nodal diameters, but the method shows large discrepancies with Method 3 (which was shown to be exact). The discrepancies with the CNCMS can therefore be explained by the particular excitation used. The simulations illustrated in Figure 7 are obtained with an excitation whose amplitude is large and whose wave shape corresponds to a low wave number h ex = 3 for which the motion of the disk is important. The authors chose on purpose this set of parameters to show that the new method is more general than the CNCMS method. However note that when the condition for which the CNCMS works well (low nonlinearity and/or high excitation wave number), this method is very efficient and can be applied for both tuned and mistuned structures (as shown in [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF] and [START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF] for a realistic FEM). Through the exactness of the method seen in Figure 6 and the large computing gain given in Tables 5, the advantages of the new methodologies are clear. In the next simulations, only the reference solution and the solution of Method 3 will be computed.

Frictional effects modeled with the Dynamic Lagrangian Frequency Time (DLFT) method

In this Section, the same structure undergoing the same external excitation is used but the frictional effects are modeled with the DLFT procedure [START_REF] Nacivet | A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF]. This methodology avoids to regularize the Coulomb's law as previously done. It is highly efficient and largely used in the community [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF][START_REF] Krack | Vibration Prediction of Bladed Disks Coupled by Friction Joints[END_REF]. The purpose of this section is to validate the methodology for a discontinuous contact interface. Details on the DLFT algorithm are provided in Appendix B.

Figure 8 illustrates the frequency response function of the system with this configuration. Similarly to Figure 6, it shows that the new methodology perfectly matches the reference solution. This validates the analytical derivations of Section 3.1 for a non-regular frequency representation of the Coulomb's law.

Comparing the results of Figures 6 and8 gives some insights on the influence of the frictional effects modeling as only the expression of the nonlinear force has changed between these two sets of simulation.

For a high level of excitation, the differences are small since the system is in the "slip" situation. However for low values of amplitude, the system is partly stuck and the results differ because the DLFT is able to retrieve an exact case of a stuck system whereas the hyperbolic tangent function always allows some partial slip (in Figure 2 the slope of f nl at 0 is not perfectly vertical). This difference is mostly visible for 1 N for the ninth nodal diameter (its amplitude is twice higher in Figure 8 due to more coupling).

Complex external excitation forces

The purpose of this numerical example is to show the generality of Method 1 and to validate the fact that it can handle complex excitation forces (created with a combination of traveling and standing waves) as underlined in Section 3.2.

Let consider an external force composed of two standing wave excitations (h ex = 2 and h ex = 6) and a traveling wave excitation (h ex = 10). All forces have an amplitude of 5 N. The shape of this particular excitation is represented in Figure 9. For this excitation, only the 2, 6 and 10 nodal diameters must be computed (see Section 3.1 and Table 4) and the nonlinear forces must only be evaluated for 6 consecutive 445 sectors. Although such excitation was not considered in [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF], its methodology can be applied but requires to consider a large "supersector" composed of 12 sectors. Overall the new method enables a reduction of 91.7% for the size of the reference sector and a reduction of 50% for the number of unknowns and the number of sectors for which the nonlinear forces must be evaluated compared to [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF]. Figure 10 compares the results obtained with the new method and the reference solution. The accuracy is excellent and proves

Conclusion

Several reductions methods were proposed in this paper and offer large computational time saving by reducing both the model to consider and the number of equations to solve. They are all based on an analytical results which determine the coupling between nodal diameters for a regularized Coulomb's law. The difference in the methodologies lies in the level of approximation made: Method 1 is purely analytical and is thus the most general, Method 2 only requires the assumption of [START_REF] Petrov | A Method for Use of Cyclic Symmetry Properties in Analysis of Nonlinear Multiharmonic Vibrations of Bladed Disks[END_REF] and finally Method 3 requires a specific solution shape. These methodologies have been validated numerically with a regularized Coulomb's law, but also with the DLFT algorithm which is a non-regular frequency representation of the Coulomb's law. Overall, the new methods have shown really interesting applicability for industrial problems, especially Method 1. Different assumptions usually made in other reduction approaches may lack dynamic representativity and are valid only for a particular range of excitation frequencies. The proposed methodologies seem more robust because they are based on exact analytical development.

Future work would consist in extending this reduction methodology to mistuned structures. These systems are currently widely studied in the scientific community and present a numerical challenge as the cyclic symmetry property no longer holds. to obtain the interacting diameters. where ε is a penalty coefficient (the choice of this parameter is provided in [START_REF] Charleux | Etude des effets de la friction en pied d'aube sur la dynamique des roues aubagées[END_REF]). ṽr represents the harmonics of the spectral relative velocity of the system's interface and must satisfy Coulomb's law. The prediction is projected in the physical and time domains (see Equation ( 55) and the AFT procedure [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF]).

Two cases must be considered: if f nl,u < µ f nl,N (sticking case), then the prediction is correct (there is no relative velocity) and hence

f nl = f nl,u (B.5)
otherwise if f nl,u ≥ µ f nl,N (sliping case), then the prediction is incorrect and must be corrected by

f nl,x = f nl,u 1 - µ f nl,N f nl,u . (B.6)
Finally the nonlinear forces are equal to

f nl = f nl,u -f nl,x . (B.7)

Figure 1 :

 1 Figure 1: Representation of a cyclic structure with N = 8 sectors.

Figure 2 :

 2 Figure 2: Representation of friction law; ( ): exact Coulomb's law; ( ): law regularized with an hyperbolic tangent function.

  Example of a standing wave distribution for a time period. Example of a standing wave distribution at a given time (t = 0).

Figure 3 :

 3 Figure 3: Representation of a standing wave excitation h ex = 3, N = 24.

2 N 1 = 7 ( 1 N 2 N 1 1 N 2 N 1 = 8 (

 2171211218 a) Case 1 with N 1 odd. = 10 (b) Case 3 with N 1 even. c) Case 2 with N 1 even and divisible by 4.

Figure 4 :

 4 Figure 4: Numbers of sectors to consider for different N 1 -values in the case of a standing wave excitation. 1 st line: spatial distribution of the excitation over N 1 sectors. 2 nd line: cyclic representation of these distributions.

3 .

 3 The external force is applied to the tip mass m 1 of each sector. The full structure is composed of 24 identical sectors. The following Figures, results are directly represented for the spectral components and only non-zero components are plotted.

Figure 5 :

 5 Figure 5: Representation of the test case.

Figure 6 :

 6 Figure 6: Frequency response function under a standing wave excitation of various amplitudes (1 N to 5 N). Frictional effects are modeled with a regularized Coulomb's law. ( ): reference solution; ( ): Method 1; ( ): Method 2; ( ): Method 3; ( ): Petrov's method.

Figure 7 :

 7 Figure 7: Frequency response function for a 5 N standing wave excitation. Frictional effects are modeled with a regularized Coulomb's law. ( ): Method 3; ( ): CNCMS procedure.

Figure 8 :

 8 Figure 8: Frequency response function under a standing wave excitation of various amplitudes (1 N to 5 N). The DLFT procedure is used to model frictional effects. The legend matches the one in Figure 6.

  using a Schur condensation, as explained in[START_REF] Poudou | Hybrid Frequency-Time Domain Methods for the Analysis of Complex Structural Systems with Dry Friction Damping[END_REF], givesZr cr + c f nl = c f r,ext , (B.1)where Zr is the relative flexural rigidity matrix of the system, cr denotes the relative harmonics components of the spectral displacement, c f nl corresponds to the harmonics of the nonlinear forces projected on the spectral basis, and c f r,ext contains the harmonics of the reduced spectral external forces. In the case of a frictional contact, the nonlinear forces are sought as c f nl = c f r,ext -Zr cr + ε (∇c r -ṽr ) , (B.2)

  The nonlinear force is separated into two parts: the prediction of the nonlinear forces c f nl,u = c f r,ext -Zr cr + ε∇c r , (B.3) and the correction of the nonlinear forces to satisfy Coulomb's law: c f nl ,x = -ε ṽr . (B.4)

Table 2 :

 2 Comparison of the different methods for a traveling wave excitation.

	Methodology	Size of the reference sector	Nonlinear solver (HBM) Number of unknowns AFT procedure	Rank
	Reference

Table 4 :

 4 Comparison of the values of N 1 and N 2 for different values of N and h ex .

	355
	12 sectors. With Method 3, only one single sector needs to be created and the nonlinear evaluation is done
	for 3 sectors.
	N			24			27
	h ex	2	3	4	5	2	3	4	5
	N 1	12	8	6	24	27	9	27	27
	N 2	3	2	2	6	14	5	14	14
	q interact (2, 6, 10) (3, 9) (4, 12) (1, 3, 5, 7, 9, 11)	0, 13	(0, 3, 6, 9, 12)	0, 13	0, 13
	N diam	6	4	4	12	27	9	27	27

Table 5 :

 5 Performance of the methods. The computation time is given for the 5 N amplitude. The simulations were run on a Intel(R)

	Core(TM) i7-7700 @ 3.6 GHz computer.

Acknowledgements

The authors would like to acknowledge the financial support of Safran Helicopter Engines.

that Method 1 can perfectly handle complex excitation and large coupling. The DLFT procedure is used to model the frictional effects. The legend matches the one in Figure 6.

Appendix

A. General algorithm to determine nodal diameter interactions

In Section 3.1 analytical derivations were proposed to obtain the interacting nodal diameters for an excitation with a specific wave number h ex . This section proposes a general algorithm enabling to determine the coupling when multiple excitations with different wave numbers are applied to the structure. Let consider P wave numbers noted q = (q i ) i∈ 1,P . Then Equation ( 27) is split into multiple sums with the different nodal diameters and becomes

with the condition of P ∑ i=1 m i = 2n -1 (m i can be even or odd). This equation can be written as

with the condition of

δ ki being odd. The space S P denotes all the possibilities for (m i ) i∈ 1,P . For instance if P = 5, then we have three sets of possibilities : all m i are odd, 3 out of the 5 m i are odd or only 1 value out of 5 is odd. It can be easily shown than the number of possibilities is equal to P 2 if P is even and P+1 2 is P is odd.

480

For the case of a single value h ex , we can determine the periodicity (noted M i ) of each nodal diameter q i through Equation ( 32). Therefore we can compute all the possibilities of condition (A.2) by taking

To facilitate the general understanding, Algorithm A.1 summarizes the different steps to follow in order Algorithm A.1: Determination of nodal diameters coupling 1: Provide N and q 2: diam ← q {Total number of nodal diameter.} 3: for i = 1 : 1 : P do 4:

Compute M i {see Equation (32)} additional combinaison← nchoosek(diam, 2k -1){We extend list combinaison with δ ki q i ; see second term of Equation (A.2) where nchoosek denotes the binomial coefficient.} 17: end for{Keep only the nodal diameter 0 ≤ q ≤ K} 18: q update ← q update ≤ K 19: Update of diam {Retrieve only new nodal diameters.} 20: q new ← setdiff(diam, q new )

B. Dynamic Lagrangian Frequency Time method

In this Appendix, we recall briefly the procedure of the DLFT method for cyclic components (its adaptation to a full physical system is similar). In the following, we will suppose that the solids are always in contact and thus only the tangential movement is considered. Employing the HBM procedure on (40), and