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

Benjamin Chouvion

Fabrice Thouverez

Nonlinear cyclic reduction for the analysis of mistuned cyclic systems

Keywords: Mistuned Structure, Friction nonlinearities, Complex nonlinear normal modes, Nonlinear reduction, Harmonic Balance Method

Predicting the vibratory response of bladed-disks in turbomachinery is of the utmost importance to design a reliable and optimized engine. Yet, such simulations are challenging mainly due to the large size of the finite-element model used to describe the system, the existence of random mistuning coming from manufacture tolerances, and the nonlinear effects arising from the different components and their coupling. As a consequence, very few reduced-order models handling nonlinear mistuned systems have yet been proposed and the existing ones try to find the best compromise between a computationally efficient simulation and a correct description of the nonlinear phenomena. In this paper, a novel approach is proposed to tackle this challenge. Its key ideas rely on the substructuring concept and the cyclic properties of the structure. A reduction basis composed of cyclic complex nonlinear normal modes is created on each sector of the system leading to a compact nonlinear superelement per sector. The system is then assembled and a synthesis simulation is performed. Due to the main concepts employed, the method is referred to as the Substructuring method based on Nonlinear Cyclic Reduction (SNCR). It is validated on a finite element model of a tuned, a randomly mistuned and an intentionally mistuned bladed-disks. It will be shown to be fast and accurate despite strong nonlinearities. Due to its flexibility and efficiency, the method is expected to have a wide range of possible applications within the community.

Introduction

Cyclic structures are widely encountered in engineering applications: wind turbines, power stations, turboengines, etc. Determining the forced response of such systems using their cyclic symmetry properties [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF] is trivial if the structure is linear: the system of equations of motion is projected on the Fourier basis and gets decoupled on its different nodal diameters [START_REF]NASTRAN cyclic symmetry capability[END_REF][START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF]. However, when nonlinearities are taken into account (such as friction or geometrical effects), the system of equations obtained after projection gets once again coupled and may be too large to be easily solved. Many scientists have studied the combination of the mathematical property of cyclic structures and nonlinearities [START_REF] Vakakis | Dynamics of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Georgiades | Modal analysis of a nonlinear periodic structure with cyclic symmetry[END_REF][START_REF] Grolet | Free and forced vibration analysis of a nonlinear system with cyclic symmetry: Application to a simplified model[END_REF]. Based on different assumptions, some numerical strategies have been proposed to compute efficiently the nonlinear response of cyclic systems [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF][START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF].

However, the cyclic symmetry property is generally not applicable to real-life turboengines because they contain small manufacturing imperfections. In turbomachinery, these imperfections that are tolerated in the design stage are known as random mistuning. It has been shown in linear systems that the largest vibration amplitude of mistuned structures may vary significantly with a small fluctuation of mistuning and is, for most cases, higher than the one of its tuned counterpart [START_REF] Castanier | Modeling and analysis of mistuned bladed disk vibration : Status and emerging directions[END_REF]. Furthermore, mistuning may give rise to localization phenomena [START_REF] Wei | Localization phenomena in mistuned assemblies with cyclic symmetry part I: free vibrations[END_REF] in which the energy of the system is mostly concentrated within a small part of the structure.

In [START_REF] Castanier | Using intentional Mistuning in the design of turbomachinery rotors[END_REF], it was shown that introducing intentional mistuning allows better control of the system dynamics.

Strategies to account for mistuning have been developed and many of them employ the cyclic properties of the structures in their reduction basis [START_REF] Yang | A reduced-order model of mistuning using a subset of nominal system modes[END_REF][START_REF] Lim | Compact, generalized component mode mistuning representation for modeling bladed disk vibration[END_REF][START_REF] Mbaye | A reduced-order model of detuned cyclic dynamical systems with geometric modifications using a basis of cyclic modes[END_REF][START_REF] Madden | Reduced-order modeling approach for blisks with large mass, stiffness, and geometric mistuning[END_REF]. These reduced-order methodologies (ROM) allow fast-computation of the system dynamics and make possible detailed parametric studies.

Combining both mistuning and nonlinear effects is the new challenging problem [START_REF] Mitra | Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces[END_REF] whose main difficulty is to create a ROM that provides an accurate description of the system while being sufficiently reduced to allow parametric studies. Szpempliska proposed a method to synthesize a frequency forced response from a single nonlinear normal mode (NNM) [START_REF] Szemplińska-Stupnicka | The modified single mode method in the investigations of the resonant vibrations of non-linear systems[END_REF]. The NNM concept was first defined by Rosenberg [START_REF] Rosenberg | The normal modes of nonlinear n-degree-of-freedom systems[END_REF] to describe the property of nonlinear systems and to gain insight into their behavior. Different definitions of the NNM have since been proposed taking into account different mechanical features such as damping [START_REF] Shaw | Normal Modes for Non-Linear Vibratory Systems[END_REF] and friction nonlinearities [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF][START_REF] Krack | Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept[END_REF]. A review of the computation and applicability of these NNMs was proposed by Kerschen et al. [START_REF] Kerschen | Nonlinear normal modes, Part I: A useful framework for the structural dynamicist[END_REF][START_REF] Peeters | Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques[END_REF]. In addition to the work of Swpempliska [START_REF] Szemplińska-Stupnicka | The modified single mode method in the investigations of the resonant vibrations of non-linear systems[END_REF], nonlinear forced response synthesis has also been used, for instance in [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF]. Based on this method, Joannin et al. proposed the Component Nonlinear Complex Mode Synthesis method [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF][START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF]. This also relies on the use of NNMs but those are calculated with fixed cyclic boundaries. The modes are exploited to create a reduction basis and to simulate the nonlinear forces by substitution. Unfortunately, it may lack accuracy when the nonlinearities strongly impact the cyclic boundaries (see 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]). Other strategies have been proposed recently [START_REF] Tang | Reduced-order models for blisks with small and large mistuning and friction dampers[END_REF][START_REF] Pourkiaee | A reduced order model for nonlinear dynamics of mistuned bladed disks with shroud friction contacts[END_REF][START_REF] Mashayekhi | Hybrid reduction of mistuned bladed disks for nonlinear forced response analysis with dry friction[END_REF], but these require to evaluate the nonlinear force for all sectors of the system with an Alternating Frequency Time procedure [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF]. This operation is numerically expensive and the methodologies seem applicable neither if one intends to perform parametric studies (such as a Monte-Carlo simulation [START_REF] Mitra | Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces[END_REF]) nor if the contact interface contains many degrees of freedom (DOFs).

The purpose of the work presented in this paper is to propose a new efficient method to tackle nonlinear mistuned cyclic systems. The manuscript is constructed as follows: Section 2 shows some limitations of current strategies that can be used to model mistuned cyclic structures, and illustrates the need for a new ROM. Section 3 presents the key concept of the new approach: a cyclic synthesis based on the cyclic complex NNMs of the structure. Section 4 then fully details the new method named "Substructuring method based on Nonlinear Cyclic Reduction" (SNCR). Finally, Section 5 provides extensive numerical validation on a finiteelement model. In this paper, a bladed-disk model is used as a test case but the proposed methodology could be employed for any cyclic structures such as turbines, pumps, etc. Moreover, this paper handles contact interfaces with friction and possible separation but the methodology could also be extended to account for geometrical nonlinearities.

Strategies to handle nonlinear mistuned structures: application and limitations

This section presents two existing strategies that will be employed later in this paper to benchmark the new methodology developed in Section 4. A reference methodology is first proposed and then a ROM, called the CNCMS [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF], is presented. The CNCMS is chosen here because it shares similar features with the new methodology presented in Section 4, such as the nonlinear modal synthesis. The purpose of this section is also to show some limitations of these possible strategies and the need to develop a more robust and less time-consuming method.

Common step: linear Craig-Bampton component mode synthesis

Turboengineers are nowadays designing intentionally mistuned (large level of mistuning) engines in order to overcome the difficulties associated with random mistuning and especially the energy localization on specific sectors of the cyclic system [START_REF] Wei | Localization phenomena in mistuned assemblies with cyclic symmetry part I: free vibrations[END_REF][START_REF] Castanier | Using intentional Mistuning in the design of turbomachinery rotors[END_REF]. An intentionally mistuned system is composed of M different kinds of sector (noted (S i ) i∈ 1,M ). Each of these sectors S i occurs m i times in the structure following a

M ∑ i=1 m i = N, (1)
where N is the total number of sectors (an example of intentional mistuning pattern is provided in Section 5.2.2). A sector can be composed of different parts such as a blade, a portion of a disk, underplatform dampers, shrouds,... Components Mode Synthesis (CMS) methods [START_REF] Tran | Component mode synthesis methods using partial interface modes: Application to tuned and mistuned structures with cyclic symmetry[END_REF] are widely used to reduce the size of a model while retaining sufficient information for its accurate description. A Craig-Bampton (CB) CMS [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] is here performed on the "blade" and the "portion of disk" components forming each kind of sector A standard CB reduction matrix is obtained and applied to both the blade and the portion of the disk. The resulting reduced mass M and stiffness matrices K of this superelement have a specific layout [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. Associated with the blade, they can be written as:

(S i ) i∈ 1,
M blade =      M nl,nl M nl,tip M nl,η M tip,nl M tip,tip M tip,η M η,nl M η,tip I      K blade =      K nl,nl K nl,tip 0 K tip,nl K tip,tip 0 0 0 Ω 2 η      , (2)
where Ω η denotes the internal natural frequencies matrix derived from the CB procedure. The matrix I represents the identity matrix while 0 is the zero matrix. The layout is similar for the disk portion, simply replace the subscript "tip" with "b" for the boundary DOFs. This first CB procedure is performed on each different kind of sector S i defining the intentional mistuning pattern. This reduction step (illustrated between States 0 and 1 in Figure 1) is common to all methodologies presented afterward.

In this paper, random mistuning is introduced by modifying the internal vibration frequencies Ω η (appearing in the reduced stiffness matrix, see Equation ( 2)) in this initial CB procedure. This step is also common to all methodologies and defines State 1b in Figure 1.

Reference methodology

For the associated linear system (stuck case), the different sectors are assembled along their boundaries as illustrated on the left column in Figure 1, and no further reduction is applied. This methodology is referred

+ ∀i ∈ 1, M • + ∀i ∈ 1, M • + ... + • + ... + • State 0 State 1

State 1b

Craig-Bampton

Blades random mistuning

Assembling the disk CB on the disk to as the full model.

For a full realistic 3D model, such as the one presented in Section 2.4, using this method without further reduction is not applicable in a nonlinear analysis. A second CB reduction is therefore performed on the disk to remove the boundary DOFs. Although this reduction step creates some additional approximations, this procedure will still be referred to as the reference solution. It is illustrated in the middle column of Figure 1.

Finally, the nonlinear system to solve becomes

M 360 ü360 + C 360 u360 + K 360 u 360 + f nl,360 = f ext,360 , (3)
where M 360 is the mass matrix of the full system after the two stages of CB reduction (and similarly for the stiffness and damping matrices). A modal damping of 5 × 10 -4 is applied to the structure. The vectors u 360 , f nl,360 and f ext,360 represent respectively the displacement, the nonlinear forces and the external forces for the entire system. The external forces are supposed periodic in time. The force f j,ext applied on a particular sector j can be localized, or can follow a specific wave pattern and be defined with respect to the force applied on the first sector. This excitation pattern can then be either a traveling (forward or backward) wave pattern

f j,ext = f 1,ext t ± 2πh ex ( j -1) ωN , j ∈ 1, N , (4)
or a standing wave pattern

f j,ext = 1 2 f 1,ext t - 2πh ex ( j -1) ωN + f 1,ext t + 2πh ex ( j -1) ωN , j ∈ 1, N , (5)
where h ex is the excitation wave number, ω is the excitation frequency that is usually the relative speed between two stages of the engines, and j is the number of sectors.

Different methodologies can be employed to solve (3), as proposed in [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. In the case of periodic external forces, the Harmonic Balance Method (HBM) is commonly used by the community [START_REF] Krack | Theory of harmonic balance[END_REF] as a solution procedure. In the HBM, the solution of (3) is sought after as

u 360 = 1 2   c 360,0 + N h ∑ n=1 c 360,n e inωt   + c.c, (6)
where (c 360,n ) n∈ 0,N h are the harmonic coefficients of the displacement and c.c denotes the complex conjugate terms. N h is the number of harmonics retained in the approximation. Substituting u 360 in (3) by ( 6) and projecting the system on the exponential basis e niωt n∈ 0,N h with the following scalar product,

f , g = ω π 2π ω 0 f (t) ḡ (t) dt, ( 7
)
where ḡ is the complex conjugate of g, leads to the following system

Z n c n + c f nl ,n = c f ext ,n ∀n ∈ 0, N h . (8)
The dynamic rigidity matrix Z n associated with the harmonic n is equal to (niω

) 2 M 360 + (niω) C 360 + K 360 .
The projection of the nonlinear force on the exponential basis is denoted by c f nl ,n . Respectively, c f ext ,n refers to the excitation force. The nonlinear forces are evaluated using an AFT procedure [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF]. In this paper, a

Schur condensation [START_REF] Poudou | Hybrid frequency-time domain methods for the analysis of complex structural systems with dry friction damping[END_REF] is used to describe the nonlinear phenomena (frictional contact and/or separation at interfaces) in the Dynamic Lagrangian Frequency Time procedure [START_REF] Nacivet | A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF]. In the following applications, the number of unknowns thus reaches N × N r,nl × (1 + 2N h ), where N r,nl is the relative number of nonlinear DOFs of the FE model. Three harmonics are usually sufficient to describe accurately such systems without requiring too many computational resources [START_REF] Pesaresi | Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications[END_REF].

Component Nonlinear Complex Mode Synthesis (CNCMS) methodology

The CNCMS method developed by Joannin et al. [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF][START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF] was mainly based on the work of Apiwattanalunggarn et al. [START_REF] Apiwattanalunggarn | Component mode synthesis using nonlinear normal modes[END_REF] which assembles nonlinear superelements obtained from the invariant manifold approach, and the work of Krack et al. [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF] which synthesizes a frequency forced response simulation from the computation of complex nonlinear normal modes [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF]. The key concept of the method in [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF] is the substitution of the nonlinear forces directly from the complex nonlinear normal modes. By this means, it avoids employing a time-consuming AFT procedure during the synthesis process.

The CNCMS combines these different methodologies. It is composed of three reduction steps. The first step corresponds to the linear CB reduction presented in Section 2.1. Then, for each different sector

(S i ) i∈ 1,M ,
a nonlinear CB reduction is performed. It mainly consists in evaluating a complex nonlinear normal mode [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] with clamped boundaries and computing static modeshapes by setting successively the displacement of each boundary DOFs to unity. The unknowns of the model after this second reduction are the boundary DOFs and the control coordinates of each sector. The number of boundary DOFs is usually large and thus a third and final reduction is employed: the interface modes [START_REF] Tran | Component mode synthesis methods using partial interface modes: Application to tuned and mistuned structures with cyclic symmetry[END_REF]. This step may however be costly as the matrices associated with the boundary can be large. Finally, the unknowns of the system are the control coordinates (q j ) j∈ 1,N , which control the complex NNMs, and the interface modes control coordinates p. The equations of motion obtained are solved by the HBM as explained in Section 2.2. The first harmonic coefficients of q j and p are denoted by c q j ,1 and c p,1 . Using a single harmonic approximation, the system to solve reduces to the algebraic system:

Z r, j   c q j ,1 c p,1   + c f nl,r, j = c f ext,r, j , (9)
where c f nl,r, j is the first harmonic coefficient of the reduced nonlinear forces, and similarly for the external forces. Z r, j is the reduced dynamic matrix and depends on the control coordinate q j . In practice the term c f nl,r, j is calculated by substitution via the complex nonlinear normal mode evaluated during the second stage of reduction [START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF][START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF]. Equation ( 9) is then assembled along the interfaces control coordinates to account for all sectors. A semi-analytical jacobian matrix is computed, as proposed in [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF], and allows fast iterations of the nonlinear solver.

Description of the test case

All numerical results presented in this paper are obtained from the fundamental bladed-disk sector represented in Figure 2. The model was made with the FE-software Ansys, and its fundamental sector is described by 25920 DOFs. 180 DOFs are nonlinear (contact interfaces with friction and possible separation) and located on the dovetail attachment, and 2130 DOFs are located on the cyclic boundaries. A rotational speed of 1000 rad s -1 is applied and leads to static preload, centrifugal softening and geometric stiffening [START_REF] Geradin | Mechanical vibrations: theory and application to structural dynamics[END_REF].

The gyroscopic effects are negligible for this model and are thus not accounted for. The new reduced-order methodology, presented in Section 4, could handle such effects if needed. The friction coefficient is taken equal to 0.3.

The excitation frequency range focuses on the first flexural mode. We have ensured that sufficient information were kept in the reduction steps (linear CB, complex NNM, interfaces modes) of the methodologies to accurately account for this mode. The number of internal and interface modes used are provided in After these reduction steps and a Schur condensation [START_REF] Poudou | Hybrid frequency-time domain methods for the analysis of complex structural systems with dry friction damping[END_REF], the reference model contains a total of 24 ×

Numerical simulation

An excitation of 0.25 N with a wave number equal to h ex = 3 is applied at the tip of the blades. The excitation follows either a traveling wave shape (see Equation ( 4)) or a standing wave shape (see Equation ( 5)).

No mistuning is considered at this point. The nonlinear solution is calculated with the reference method-145 ology (see Section 2.2) and the CNCMS method (Section 2.3). The sectors are expected to undergo the same displacement shifted in time [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF] when excited by a traveling wave, whereas a standing wave excitation provides different levels of energy for each sector and implies different displacement amplitudes between the blades. Figure 3 illustrates the tip amplitude of all the blades. For the linear underlying system (stuck case), the full model solution, the linear reference solution, and the CNCMS method give the same results. Only the linear response obtained from the reference method is illustrated in Figure 3. As expected, the traveling wave excitation (see Figure 3a) yields an amplitude equal for all sectors. However, the response from the standing wave excitation, given in Figure 3b, shows three different levels (one of them is negligible, below 10 -6 m). For the two highest levels, the CNCMS method is not able to describe accurately the vibration amplitude and one can observe a decrease of approximately 8% at the resonant peak.

These discrepancies come from the lack of representativeness of the nonlinear CB reduction procedure in the CNCMS method. The nonlinear basis contains modes calculated with fixed boundaries and uses only linear static modeshapes for static correction. Such a basis is expected to provide good accuracy when either the motion of the disk is small (the system mainly undergoes a blade motion mode) or when the level of nonlinearity is low and does not impact the boundary DOFs. In the case under study, the second reason is mostly responsible for the discrepancies observed.

Table 2 provides the computation time of the nonlinear simulations performed with the CNCMS and the reference methodology. In this example, the CNCMS is approximately seven times faster than the reference.

Moreover, if one is interested in performing a Monte Carlo analysis to study the impact of random mistuning (requiring thousands of simulation), the CNCMS re-uses the pre-computed NNMs and interface modes and thus only need to relaunch the synthesis procedure.

The main objective of the paper is therefore to propose a new approach that can be more robust than the CNCMS method while still offering fast simulations.

CNCMS Reference

NNM Interface Modes Synthesis

Number of unknowns 630 848 15120

Computation time 10 min 11 min 40 min 7 h

Table 2: Computation times for the nonlinear simulations of the system under a traveling wave excitation. The simulations were run on a Intel(R) Core(TM) i7-7700 @ 3.6 GHz computer.

Nonlinear cyclic modal synthesis 170

The new methodology developed in this paper employs both the concept of substructuring, as explained in Section 2.1, and the concept of cyclically symmetric complex nonlinear normal modes. This latter key idea is first explained and used in a nonlinear cyclic modal synthesis to simulate the response of tuned cyclic structures.

Cyclic complex nonlinear normal mode

One property of cyclic symmetric structures, i.e. without mistuning, is that the displacement of its j th sector can be written under the Fourier form [START_REF] Valid | Théorie et calcul statique et dynamique des structures à symétries cycliques[END_REF]

u j = N-1 ∑ k=0 ũk e -i( j-1)kα , j ∈ 1, N , (10) with ũN
-k = ūk , k ∈ 1, N -1 (11)
and α = 2π N . ũk are the spectral components and knowing them for k ∈ 0, K (where

K = N 2 if N is even and K = N-1
2 otherwise) is sufficient to recover the physical displacement u j . In practice, the components ũk are decomposed into the internal DOFs (subscript i), and the left and right boundary DOFs (subscripts L and R). It respects the following property [START_REF]NASTRAN cyclic symmetry capability[END_REF][START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF]:

ũk =      ũk,i ũk,b,L ũk,b,R      =      I 0 0 I 0 e -ikα I      Bk   ũk,i ũk,b,L   , (12)
where the left and right boundaries are linked by a phase change function of the nodal diameter k. By means of the matrix Bk k∈ 0,K and a projection on the Fourier basis (Equation ( 10)), we can rewrite the system of equations of motion for both linear [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF] and nonlinear problems [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF][START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]. This leads to fully decoupled equations for linear systems. However, for nonlinear systems, the different nodal diameters interact (the interacting diameters define the interaction space T interact [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF][START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]) and the system of equations becomes,

Mk ük + Ck uk + Kk ũk + fnl,k = fext,k , k ∈ T interact , (13)
where Mk , Ck and Kk denote the mass, stiffness and damping matrices of the k-th nodal diameter. The vectors ũk , fnl,k and fext,k represent the displacement, the nonlinear forces and the external forces on the k-th nodal diameter. T interact can be taken by default as the whole space 0, K but it might be possible to reduce it depending on the situation [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF][START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]. Solving [START_REF] Lim | Compact, generalized component mode mistuning representation for modeling bladed disk vibration[END_REF] with the term fext,k allows to compute the frequency forced response of the system. The solution of the associated autonomous system (setting fext,k = 0 in (13)) gives 180 its cyclic nonlinear normal modes.

In this Section, we first focus on the autonomous system for a specific nodal diameter (k = h ex ) in order to study one of its modes. In linear systems, the eigenfrequency of a degenerated nodal diameter k / ∈ 0, N 2 is double and can be described by a mode pair that consists of a forward and a backward traveling modes [START_REF] Mitra | Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces[END_REF].

Based on this linear description, we propose here to define the cyclic nonlinear normal modes as a sum of a forward and backward traveling modes. The computation of a forward traveling cyclic nonlinear normal modes for a particular nodal diameter h ex is explained next. The computation of a backward traveling mode can be obtained similarly. Equation ( 13) is solved with fext,k equal to 0 (the static preload is considered but the dynamic excitation is set to 0) and by assuming that the solution can be written as a forward traveling wave [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF]

u j (t) = u 1 t - α ( j -1) h ex ω . ( 14
)
Different methodologies have been proposed to compute the nonlinear modes of a system undergoing friction nonlinearities. In [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF], the dissipation introduced by friction was characterized by a parameter β , whereas it was described in [START_REF] Krack | Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept[END_REF] with an additional damping matrix. As one is here interested in synthesizing a frequency forced response based on complex nonlinear modes, it seems better-suited to parametrize the 185 friction by a single parameter β instead of a full matrix. The definition of [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF] is then adopted here. Using this definition, the solution on the first sector is sought after as a modified HBM form:

u 1 (t) = 1 2   c 0 + N h ∑ n=1 c n e n(-β +iω)t   + c.c. ( 15
)
The displacement of the remaining sectors is evaluated with Equation ( 14). By comparison with the traditional HBM (see Equation ( 6)), one can see that this expression adds an exponential term (e -βt ) related to the dissipation effects. In the case of friction nonlinearities, the characteristic time associated with β is much longer than the period associated with ω. Therefore, it may be assumed, as in [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF], that β ω. Combining the assumption of a forward traveling wave shape (Equation ( 14)) with the solution (15) leads to a correspondence between the n-th harmonic and the k-th nodal diameter (governed by the number of sectors and the nodal diameter under study h ex ) [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF][START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]. This relationship imposes the displacement of the nodal diameter k to be governed by several harmonics (whose space will be noted N k ), whereas a given harmonic n is only present in a single nodal diameter expansion. The spectral components are hence equal to

ũk =      √ N 2 cn e nλ - hex t or √ N 2 c n e nλ + hex t if k ∈ 1, M , n ∈ N k √ N 2 c n e nλ + hex t + cn e nλ - hex t if k ∈ 0, N 2 , n ∈ N k , (16)
where λ ± h ex = (-β ± iω) and M = N 2 -1 if N is even and M = N-1 2 otherwise. The sign of the scalar λ h ex depends on the harmonic n which may lead to a nodal diameter aliasing.

As in the traditional HBM, the final system to solve is obtained by substituting ( 16) in (13) (with fext,k = 0) and projecting the system on the exponential basis e niωt n∈ 0,N h with the scalar product defined in [START_REF] Petrov | A method for use of cyclic symmetry properties in analysis of nonlinear multiharmonic vibrations of bladed disks[END_REF].

Since it is assumed that β ω, the projection of e -nβt on the exponential basis e niωt n∈ 0,N h is assumed to be equal to 1, and we finally have the following equations:

         λ - h ex 2 Mh ex + λ - h ex Ch ex + Kh ex √ N 2 c1 + fnl,h ex , e -iωt = 0 n = 1 nλ ± h ex 2 Mk + nλ ± h ex Ck + Kk √ N 2 cn + fnl,k , e ±inωt = 0 ∀n ∈ 0, N h and n = 1 (17a) (17b)
where fnl,k , e ±inωt are the projections of the spectral nonlinear forces on the exponential basis. Details of these derivations can be found in [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF] for a similar case. The nonlinear forces will be evaluated with the AFT procedure [START_REF] Cameron | An alternating frequency/time domain method for calculating the steadystate response of nonlinear dynamic systems[END_REF].

The unknowns of the system (17) are the harmonic coefficients (c n ) n∈ 0,N h , the frequency ω and the coefficient β . This system is underestimated and two additional real equations must be provided. Imposing the amplitude and phase of the first harmonic coefficient of the observation DOF (here taken as the tip of the blade) can be used to supplement the system of equation [START_REF] Laxalde | Complex non-linear modal analysis for mechanical systems: Application to turbomachinery bladings with friction interfaces[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF]:

ℜ c1,tip = x amp , (18a)
ℑ c1,tip = -x amp (18b)
where ℜ and ℑ represent the real and imaginary parts. The term x amp ∈ 0, X f is a parameter that will be gradually increased up to X f during the evaluation of the nonlinear normal mode. The choice of the phase equation (18b) is arbitrary. Solving simultaneously [START_REF] Szemplińska-Stupnicka | The modified single mode method in the investigations of the resonant vibrations of non-linear systems[END_REF] and [START_REF] Rosenberg | The normal modes of nonlinear n-degree-of-freedom systems[END_REF] gives the cyclic complex NNM with traveling wave shape associated with the nodal diameter h ex .

In the following applications, the pair c1 , λ - h ex represents the discretized NNM (mode shape c1 , modal frequency ω and damping coefficient β ) and only this information is retrieved from the solution of the system. For each level of x amp , c1 is normalized by c1,tip to give the normalized nonlinear mode shape Φ - h ex :

Φ - h ex = c1 c1,tip = c1 x amp (1 -i) , x amp ∈ 0, X f . (19)
Φ - h ex is numerically evaluated for discretized values of x amp (a few hundred is usually sufficient to retrieve the whole NNM behavior). To describe the motion of the system for any level of amplitude, a linear interpolation of Φ - h ex is performed. The eigenvalues of the NNM, λ - h ex , are also interpolated. In the end, Φ - h ex and λ - h ex , are no longer discretized but are function of q tip (where q tip = c1,tip = √ 2x amp is the amplitude of the tip of the blade). For better readability the dependency of Φ and λ on the tip displacement will be omitted. A similar development can be performed to obtain the backward traveling NNM defined with its mode shape Φ + h ex and eigenvalues λ + h ex .

Cyclic synthesis

The concept of cyclic complex nonlinear normal mode explained in Section 3.1 will be largely used in Section 4 to create a reduced order model of any mistuned structure. It is first used here to simulate the forced response of a cyclic tuned structure.

Let assume that the system is excited by a forward traveling wave of frequency ω and wave number h ex (h ex ≤ K). The frequency range of the excitation is supposed to be centered around a specific mode p. Based on [START_REF] Lim | Compact, generalized component mode mistuning representation for modeling bladed disk vibration[END_REF], the assumption [START_REF] Madden | Reduced-order modeling approach for blisks with large mass, stiffness, and geometric mistuning[END_REF] and the HBM, one must solve (full details can be found in [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]):

         (-iω) 2 Mh ex + (-iω) Ch ex + Kh ex √ N 2 c1 + fnl,h ex , e -iωt = fext,h ex , e -iωt n = 1 (-inω) 2 Mk + (-inω) Ck + Kk √ N 2 cn + fnl,k , e -inωt = 0 ∀n ∈ 0, N h and n = 1. (20a) (20b)
To solve Equation ( 20), we will use the strategy proposed by Krack et al. [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF] (although the method was explained for a non-cyclic system, it can be easily extended to a cyclic one). The method consists in employing the cyclic complex nonlinear normal mode, solution of Equations ( 17) and ( 18) (for the nodal diameter h ex and the mode p). The Equation (20b) does not contain an excitation term and is supposed to be verified by the solution of (17b). It then remains only Equation (20a) to solve. The term c1 is substituted by Φ - h ex q - h ex , where q - h ex is the control coordinate of the NNM. In this specific case where only one NNM is considered we have q tip = q - h ex . Equation (20a) is then pre-multiplied by Φ - h ex and becomes

Φ - h ex T (-iω) 2 Mh ex + (-iω) Ch ex + Kh ex √ N 2 Φ - h ex q - h ex + fnl,h ex , e -iωt = fext,h ex , e -iωt = 0. ( 21
)
The term fnl,h ex , e -iωt is substituted by the NNM (see Equation (17a)):

fnl,h ex , e -iωt = -λ - h ex 2 Mh ex + λ - h ex Ch ex + Kh ex √ N 2 Φ - h ex q - h ex . (22)
Finally the system to solve becomes

Φ - h ex T (-iω) 2 -λ - h ex 2 Mh ex + (-iω) -λ - h ex Ch ex √ N 2 Φ - h ex q - h ex -fext,h ex , e -iωt = 0. ( 23
)
This nonlinear system of equations has only two unknowns: the real and imaginary part of q - h ex , and no evaluation of the nonlinear forces is required. Therefore solving [START_REF] Peeters | Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques[END_REF] is relatively fast with an iterative numerical algorithm. Although the methodology was recalled for a single mode, the basis Φ may contain multiple modes in order to account for a wider range of excitation frequency. This method is expected to give very accurate results as long as there is no internal resonance within the modes of a given nodal diameter [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF][START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF]. If one is interested in a backward traveling wave excitation, the same procedure can be applied and the final equation will be similar to [START_REF] Peeters | Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques[END_REF] with the quantities Φ + h ex and λ + h ex . Let consider a practical case to illustrate the effectiveness of this nonlinear cyclic synthesis. The tip of the blades of the structure represented in Figure 2 is excited with a forward traveling force with an amplitude of 0.1, 0.15, 0.2 or 0.25 N and a wave number h ex = 3. The reference solution is computed (see Section 2.2) as well as the synthesis procedure proposed in this section. As seen and explained in Section 2.5, for a traveling wave excitation, the sectors will all undergo the same maximal amplitude shifted in time. As a consequence, the results of Figure 4 only represent the displacement of one of its sector. Figure 4 shows that the results of the cyclic synthesis are very accurate and that the backbone curve of the NNM [START_REF] Peeters | Nonlinear normal modes, Part II: Toward a practical computation using numerical continuation techniques[END_REF] approximately locates the peak of each frequency forced response. Furthermore, the cyclic synthesis procedure requires only a few seconds of simulation. The computation time is given in Table 3.

The corresponding computation of the cyclic NNM is five times longer than the fixed NNM (see Table 2)

because the boundary DOFs are kept as linear DOFs during the Schur condensation [START_REF] Poudou | Hybrid frequency-time domain methods for the analysis of complex structural systems with dry friction damping[END_REF] in the cyclic NNM 235 calculation.

Cyclic synthesis Reference

NNM Synthesis

Number of unknowns 630 2 15120

Computation time 50 min 30 s 7 h

Table 3: Computation times for the cyclic synthesis simulation with a 0.25 N amplitude excitation. The simulations were run on a Intel(R) Core(TM) i7-7700 @ 3.6 GHz computer.

The results of the synthesis procedure are very promising. However, as they are based on the assumption of a traveling wave shape, they are only valid for cyclic symmetric structure undergoing traveling wave excitations. The cyclic synthesis is unable to retrieve the different levels of amplitude of the sectors under a standing wave excitation. This is illustrated in Figure 5: the cyclic synthesis only gives a single level of energy (due to the traveling wave assumption) whereas the reference solution has got three distinct energy states.

The purpose of the next section is to present a new methodology which can handle any kind of excitations (traveling, standing, localized) but can also take efficiently into account mistuned structures. The SNCR employs as first step the linear substructuring detailed in Section 2.1. The other reduction steps are explained below and summarized in Figure 6. For each kind of sector (S i ) i∈ 1,M , one first assumes the entire structure to be perfectly cyclic with fundamental sector S i , and evaluates the associated cyclic matrices (State 2 in Figure 6). Both the forward and the backward traveling cyclic complex nonlinear normal modes are then computed, as explained in Section 3.1 (see State 3 in Figure 6). Since the spectral displacement can be written as the sum of the forward and backward modes over all modes [START_REF] Mitra | Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces[END_REF], we assume that the spectral component of a particular nodal diameter can be written on a nonlinear modal basis as the sum of forward and backward traveling cycling complex nonlinear modes modulated by different amplitudes:

ũk,i = N mode ∑ p=1 Φ - k,p,i q tip q - k,p e -iωt + Φ + k,p,i q tip q + k,p e iωt , i ∈ 1, M , (24)
where N mode is the truncated number of modes kept for the analysis and ũk,i represents the spectral displacement of the k-th nodal diameter for the type of sector i. The quantities q - k,p and q + k,p are the unknown generalized coordinates of the problem associated with the nodal diameter k and mode number p. Notice that for non-degenerated nodal diameters there is a single eigenvector per mode (with k = 0, all sectors move 250 in unison, and with k = N 2 the sectors have an alternating motion). Let consider a sector j with the intentional mistuning type i ∈ 1, M . Based on Equations ( 10) and ( 24), the displacement for this sector is

u j = N-1 ∑ k=0 ũk,i e -i( j-1)kα (25a) = N-1 ∑ k=0   N mode ∑ p=1 Φ - k,p,i q tip, j q - k,p, j e -iωt + Φ + k,p,i q tip, j q + k,p, j e iωt   . (25b)
The term e -i( j-1)kα in (25a) is directly included in the amplitudes q - k,p, j and q + k,p, j in Equation (25b). Notice that the formulation in (25b) is more flexible than [START_REF] Wei | Localization phenomena in mistuned assemblies with cyclic symmetry part I: free vibrations[END_REF] because each sector j is here controlled by its own generalized coordinates. The displacement u j must satisfy the following equation of motion for a sector j,

M j ü j + C j u j + K j u j + f nl, j + f s, j = f ext, j , (26)
where M j , C j , K j are the mass, damping and stiffness matrices of the sector j and account for both the intentional mistuning (sector j is of type i ∈ 1, M ) and the random mistuning. The term f s, j corresponds to the forces applied by the neighboring sectors on the sector j itself. The HBM with a single harmonic is applied on Equation (26) as presented in Section 2.2 and gives:

Z 1, j c 1, j + f nl, j , e iωt + f s, j , e iωt = f ext, j , e iωt , (27)
where Z 1, j is the dynamic matrix associated with the first harmonic of sector j and c 1, j is the first complex harmonic coefficient associated with u j . Based on Equation (25b), the first harmonic of u j , c 1, j , is written using the nonlinear modal basis

c 1, j = u j , e iωt = N-1 ∑ k=0 N mode ∑ p=1
Φ - k,p,i q tip, j q - k,p, j e -iωt + Φ + k,p,i q tip, j q + k,p, j e iωt , e iωt

= K ∑ k=0 N mode ∑ p=1
Φ - k,p,i q tip, j q - k,p, j e -iωt + Φ + k,p,i q tip, j q + k,p, j e iωt , e iωt

+ M ∑ k=1 N mode ∑ p=1 Φ- k,p,i q tip, j q- k,p, j e iωt + Φ+ k,p,i q tip, j q+ k,p, j e -iωt , e iωt = K ∑ k=0 N mode ∑ p=1 Φ- k,p,i q tip, j q- k,p, j + Φ + k,p,i q tip, j q + k,p, j , (28)
where the property [START_REF] Castanier | Using intentional Mistuning in the design of turbomachinery rotors[END_REF] has been employed to simplify the initial sum (from K + 1 to N -1) into the sum from 1 to M (with complex conjugates terms). Due to the normalization of the NNMs (see Equation ( 19)), the blade tip amplitude of sector j is expressed with

q tip, j = N-1 ∑ k=0   N mode ∑ p=1 q- k,p, j + q + k,p, j   . (29)
In the following, the eigenvector for a given nodal diameter k, a type of sector i and a mode p that combines both forward and backward traveling modes, is noted,

Φ k,p,i = Φ- k,p,i Φ + k,p,i , (30)
where the dependency on q tip, j has been dropped for better readability. The complete basis is numerically created by arranging the NNMs one after the other

Φ i = Φ 0,1,i ... Φ h ex ,1,i ...Φ K,1,i ... Φ h ex ,p,i ... Φ K,N mode ,i . (31)
The vector of unknowns follows the same decomposition:

q k,p, j = q- k,p, j q + k,p, j T (32)
and

q j = q T 0,1, j ... q T h ex ,1, j ... q T K,1, j ... q T h ex ,p, j ... q T K,N mode , j T . ( 33
)
Equation ( 28) is finally written as

c 1, j = Φ q tip, j q j . ( 34
)
To account for the impact of the neighbouring sectors on sector j itself, one must assemble all sectors. The displacement is therefore decomposed into internal DOFs (subscript i), and the boundaries DOFs (subscript b, L and b, R for the left and right boudaries),

c 1, j =      c 1, j,i c 1, j,b,L c 1, j,b,R      =      Φ i,i q tip, j Φ i,b,L q tip, j Φ i,b,R q tip, j      q j . ( 35
)
This equation is then modified to control independently the internal, left and right DOFs,

c 1, j =      c 1, j,i c 1, j,b,L c 1, j,b,R      =      Φ i,i q tip, j 0 0 0 Φ i,b,L q tip, j 0 0 0 Φ i,b,R q tip, j      Φ i,f      q j q j,b,L q j,b,R      . ( 36
)
The previous operation adds flexibility to the system with the additional unknowns q j,b,L , q j,b,R and enables the full structure assembly. Notice that the tip amplitude of sector j is still controlled by the quantity q j through Equation ( 29).c 1, j is substituted by ( 36) in ( 27) and then the system is pre-multiplied by Φ T i,f to give

Z r,j      q j q j,b,L q j,b,R      + Φ T i,f f nl, j , e iωt + Φ T i,f f s, j , e iωt = Φ T i,f f ext, j , e iωt , (37)
The reduced dynamic rigidity matrix is equal to Z r,j = Φ T i,f Z 1, j Φ i,f . When solving [START_REF] Nacivet | A dynamic Lagrangian frequency-time method for the vibration of dry-friction-damped systems[END_REF], the different f s, j terms cancel each other out when assembling each sector on the q j,b,L , q j,b,R control coordinates. As presented in Section 3.2, the nonlinear forces are substituted with the NNMs by using u j defined in Equation [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF] and with the assumption that the cyclic nonlinear forces are decoupled one from another:

f nl, j , e iωt = f nl (u j ) , e iωt = f nl N-1 ∑ k=0 N mode ∑ p=1 Φ - k,p,i q tip, j q - k,p, j + Φ + k,p,i q tip, j q + k,p, j , e iωt = N-1 ∑ k=0 N mode ∑ p=1 f nl Φ - k,p,i q tip, j q - k,p, j , e iωt + N-1 ∑ k=0 N mode ∑ p=1 f nl Φ + k,p,i q tip, j q + k,p, j , e iωt . ( 38
)
For each nodal diameter and mode, the term f nl Φ - k,p,i q tip, j q - k,p, j , e iωt is substituted by Equation (17a) (and similarly for the backward traveling NNM). This gives:

f nl, j , e iωt = - K ∑ k=0 N mode ∑ p=1 λ - k,p 2 Mk + λ - k,p Ck + Kk √ N 2 Φ- k,p,i q- k,p, j - K ∑ k=0 N mode ∑ p=1 λ + k,p 2 Mk + λ + k,p Ck + Kk √ N 2 Φ + k,p,i q + k,p, j , (39)
where the dependency of Φ and λ on q tip, j has been omitted.

In terms of numerical implementation, one creates for each mode p and nodal diameter k the vector

z - nl, j = -λ - k,p 2 Mk + λ + k,p Ck + Kh ex √ N 2 Φ k,p,i (40a)
z + nl, j = -λ + k,p 2 Mk + λ + k,p Ck + Kk √ N 2 Φ k,p,i (40b)
and the different vectors are concatenated in a matrix Z nl, j (similarly to the Equations ( 30) and ( 31)). Finally, the system to solve becomes

Z r,j      q j q j,b,L q j,b,R      + Φ T i,f Z nl, j q j = Φ T i,f f ext, j , e iωt . ( 41
)
As the quantities Z r,j , Φ T i,f and Z nl, j depends on q j (see Equation ( 29)), the system (41) must be solved with a nonlinear solver. Similarly to the CNCMS method [START_REF] Joannin | A nonlinear component mode synthesis method for the computation of steady-state vibrations in non-conservative systems[END_REF], a semi-analytical jacobian is provided to ensure fast computation of the solver.

Once the system is solved (the quantities q j , q j,b,L and q j,b,R are known for all sectors), one can recover the displacement for all DOFs. The way is to compute the amplitude at the tip of each sector (Equation ( 29)), then to evaluate the mode shapes for this value (via interpolation), and finally to employ Equation (25b) to obtain the physical displacement.

How to choose the nodal diameters to include in the nonlinear basis ?

The total number of unknowns of the system ( 41) is equal to the total number of NNMs kept in the basis multiplied by 2 (to account for the complex nature of the variables), multiplied again by 2 (for the split of unknowns introduced in Equation ( 36)), and finally multiplied by N (for the number of sectors). Knowing that the computation of a NNM is relatively long (see Table 3), one must keep the number of NNMs at the lowest to ensure fast computation. If the excitation frequency focuses on a specific family of mode (in our

• + ∀i ∈ 1, M State 1 • + ... + • + ... + • State 1b
Blades random mistuning

• ∀i ∈ 1, M Cyclic symmetry State 2 • ∀i ∈ 1, M NNM State 3 • ∀i ∈ 1, M Split of variables State 4 • • • • • • • • Assembly State 5 SNCR

22

case the first flexural mode of the different nodal diameters), the NNMs can be only computed for this family (N mode = 1 in ( 24)).

There is still to choose the nodal diameters on which one must compute the NNMs. As traveling and standing waves are the most encountered excitations in cyclic systems, the strategy adopted in this paper and advised to potential users assumes such excitations and is explained next.

1. Retrieve the wave numbers h ex of the imposed excitation forces (multiple values are possible)

2. Based on the work in [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF], compute the reduced set of nodal diameters T interact which can get coupled by this particular excitation1 .

The basis is then formed by the NNMs of these nodal diameters. Depending on the number of sectors and the wave numbers h ex , the basis may still be large. Finding a way to reduce it further is interesting as it would diminish the number of the NNMs to be calculated and reduce the final computation time of the SNCR. The following optional procedure is then proposed.

1. For each kind of sector i ∈ 1, M , compute the reference solution for the cyclic symmetric structure of type i under a standing wave excitation with the wave numbers h ex . The reference solution can either be obtained from the full system, as presented in Section 2.2, or from the methodology presented in [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF] (which was shown to be accurate for large FEM in [START_REF] Quaegebeur | Nonlinear dynamic analysis of 3d bladed-disks with frictional contact interfaces based on cyclic reduction strategies[END_REF]).

2. Create the basis Φ with the h ex values (solve [START_REF] Szemplińska-Stupnicka | The modified single mode method in the investigations of the resonant vibrations of non-linear systems[END_REF] and [START_REF] Rosenberg | The normal modes of nonlinear n-degree-of-freedom systems[END_REF] and concatenate (31))

3. Compute the frequency forced response with the SNCR using Φ as the reduction basis (solve [START_REF] Quaegebeur | Nonlinear dynamic analysis of 3d bladed-disks with frictional contact interfaces based on cyclic reduction strategies[END_REF] with the system assembled) 4. Iterate between steps 2 and 3 until convergence between the SNCR and the reference solution by supplementing the basis Φ with nodal diameters taken in T interact .

Once this basis is formed (following one of the procedures), one can launch simulations for both tuned and mistuned structures. One might argue that in the case of random mistuning, all nodal diameters get excited and thus the basis Φ must contain all NNMs. However, since the SNCR parametrizes each sector independently and that such mistuning is small, the basis is expected to be accurate enough.

Numerical validation, performance and limitations

In this section, the SNCR is validated on different cases: a cyclic symmetric structure, a randomly mistuned structure and an intentionally mistuned structure. The use of various examples demonstrates the possible wide applicability of the method. In these examples, the gyroscopic effects are negligible and are thus ignored. As a consequence, the backward and forward traveling cyclic normal nonlinear modes are complex conjugates:

Φ + k,p,i = Φ- k,p,i , (42)
and hence only the latter are computed.

Numerical validation for a cyclic symmetric structure

The cyclic symmetric structure is such as the one presented in Sections 2.5 and 3.2. The method has been validated for different values of h ex for both traveling and standing excitations; however, for brevity, only results associated with a standing wave excitation with two values of h ex (1 and 3) are here presented.

The standing wave excitation was shown to be more challenging (the sectors undergo different levels of amplitude, see Section 2.5 and more coupling occurs between the nodal diameters [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF]). There is no particular difficulty in the case of a traveling wave excitation. The excitation amplitude is taken equal to 0.25 N.

Let first consider the case h ex = 1. Based on the work in [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF] (and with a structure of 24 sectors), we can predict that coupling will occur between the nodal diameters {1, 3, 5, 7, 9, 11} and the 0 nodal diameter excited by the static loading due to an imposed rotation. The reduction basis initially used is {0, 1}, composed of the zero and first nodal diameters (the numbers inside the brackets denote the nodal diameters kept for the reduction). In this section, only the first flexural mode of the different nodal diameters is studied (in Equation [START_REF] Krack | A method for nonlinear modal analysis and synthesis: Application to harmonically forced and self-excited mechanical systems[END_REF], N mode = 1).

Figure 7 shows the amplitude of the tip of all sectors as a function of the excitation frequency, calculated for the linear and nonlinear systems, and with different methodologies. The linear response of the SNCR with basis {0, 1} gives the same response as the reference solution (as expected since only the first nodal diameter is excited and there is no coupling). This result is not provided in Figure 7a. The nonlinear response of the SNCR in Figure 7a shows some discrepancies with the nonlinear reference solution. It means that the {0, 1} basis is not satisfactory. It is then supplemented with the third nodal diameter and is now noted {0, 1, 3}. Figure 7b illustrates the corresponding results. The amplitude of all sectors is now correctly obtained and shows that convergence is reached. Similar results are found with h ex = 3. Figure 8a represents the response obtained with the basis {0, 3}.

Adding the ninth nodal diameter [START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF] to form the basis {0, 3, 9}, the SNCR becomes able to accurately describe the response of the system, as shown in Figure 8b. Very small discrepancies are observed which may result from the substitution of the nonlinear forces (see Equation [START_REF] Pesaresi | Modelling the nonlinear behaviour of an underplatform damper test rig for turbine applications[END_REF]) via the NNMs.

The SNCR provides very promising results as we were able to retrieve a nonlinear response under a standing wave excitation at a low cost. Moreover, the SNCR is more accurate than the CNCMS (a direct comparison can be made between Figure 3b and Figure 8b). Details on the performance of the method will be provided in Section 5.3.

To complete the study on the tuned structure, a localized excitation is now applied on the first sector. For this specific excitation, the strategy given in Section 4.2 is no longer applicable because all nodal diameters are coupled. The basis of the SNCR was then supplemented gradually up to {0, 1, 3, 5, 7} (the numbers were arbitrarily chosen to be odd) until the linear solution was considered to be accurate by comparison with the reference solution. Converged results are illustrated in Figure 9a. The nonlinear response of the system, using the same reduction basis, is presented in Figure 9b, and the differences between the reference solution and the ROM are shown to be small.

The next section focuses on the validation of the SNCR for mistuned structures.

Numerical validation for mistuned structures

The purpose of this section is to validate the applicability of the SNCR as a reduction method for nonlinear mistuned structures. A comparison with the CNCMS method, presented in Section 2.3, is also provided.

Random mistuning

In real engineered systems, the cyclic structure contains manufacture uncertainties [START_REF] Mitra | Dynamic modeling and projection-based reduction methods for bladed disks with nonlinear frictional and intermittent contact interfaces[END_REF] and exhibits a larger maximum displacement than its perfectly tuned counterpart (the ratio of both values is named the amplification factor, or AF). There exist different possibilities to take into account random mistuning in the model. For instance, in the CNCMS [START_REF] Joannin | Reduced-order modelling using nonlinear modes and triple nonlinear modal synthesis[END_REF], the eigenfrequency of the fixed-NNM was shifted following a normal distribution with mean value (µ) equal to the natural frequency of the linear counterpart of the NNM and a variable standard deviation (σ ) referred to as mistuning level.

In this paper, the procedure to account for random mistuning modeling consists in modifying the internal natural frequencies in the initial CB procedure (see Equation ( 2)). For each sector, each value of Ω η is chosen randomly following a normal distribution with µ taken as the natural frequency Ω η and a variable σ .

This random mistuning modeling is interesting for validation and practical applications because the random matrices are common for both the reference method and the ROM. Besides, the geometry and material of the sector are not directly modified and the computation of the structural matrices only requires the original sector matrices (offering some computation gain). In most applications, as for instance in [START_REF] Lim | Compact, generalized component mode mistuning representation for modeling bladed disk vibration[END_REF], only blade mistuning is considered and such practice is adopted in this paper. The choice of random mistuning modeling is method-dependent.

For both linear and nonlinear reference methodologies, the blades model is not further reduced after State 1 (see Figure 1) and the random mistuning matrices are introduced in the model assembly in State 4. The reference models retain thus their whole accuracy. In the case of the SNCR, see Figure 6, (respectively the CNCMS, see Figure 1), the blades are further reduced at step 2 with the calculation of the NNMs. Evaluating the NNMs for each new set of random mistuning would deteriorate the performance of the method. The basis Φ is therefore assumed not to be impacted by random mistuning. Such approximation is reasonable since the mistuning is small. The modification on the stiffness matrices of the blades is nevertheless still included in the synthesis stage (Equation ( 41) for the SNCR and ( 9) for the CNCMS). As a consequence, the differences in the results between the reference method and the SNCR may be due to two approximations:

the small impact of the mistuning on the NNMs that is not accounted for, and the reduction methodology itself (nonlinear basis Φ, substitution of the nonlinear forces, etc.).

3 is applied at the tip of the blades. The {0, 3, 9} reduction basis is taken here as suggested in Section 4.2.

Similar results are found for different values of h ex . Figure 10a and 10b represent the linear solution (stuck system) obtained with the full model (left column in Figure 1), the reference model (middle column in Figure 1) and the SNCR solution for a relatively large random mistuning pattern (equal to σ = 15%). The reference solution matches perfectly the full model 365 calculation, and it proves that the 600 internal modes kept in the nonlinear reference solution are sufficient to capture accurately the motion of the disk. For the SNCR solution, we observe that the amplitudes of the displacements are close to the reference ones although a frequency shift occurs at lower frequencies (shift of 2 Hz around 850 Hz, which corresponds to a 0.25% absolute relative error). This small discrepancy is due to the reduction procedure of the SNCR which cannot capture perfectly the flexibility of the disk. These residual flexibilities can be corrected by taking a larger base (as it was performed in Section 5.1). However, the purpose is here to show that the general behavior is correctly captured even with a very small reduced basis. It can also be noticed in Figure 10a and 10b that most of the sectors reach their maxima at higher frequencies (around 865 Hz) and the results of the SNCR are accurate in this frequency range.

Figure 10c and 10d represent the nonlinear response obtained with the different methodologies. Similarly to the linear case, a slight shift of frequency is observed around 850 Hz. However, the maximal displacement of the system (and thus the amplification factor) around 862 Hz is accurately captured. This ratio is of the utmost importance for the designers and we observe that the proposed ROM manages to evaluate it correctly.

For this specific mistuning pattern, the SNCR gives very satisfactory results. Additional simulations are provided next to fully validate the method for different σ values.

where u j,tip,ROM (respectively u j,tip,Ref ) represents the displacement of sector j at the tip of the blade evaluated with the ROM (respectively the reference solution). The error L ∞,ω is evaluated with

L ∞,ω =     max j ω j,MAX,ROM -max j ω j,MAX,Ref max j ω j,MAX,Ref     , (45)
where ω j,MAX,ROM (respectively ω j,MAX,Ref ) represents the frequency of the ROM (respectively the reference solution) for which the tip displacement of sector j reaches its maximum. In order to study the global Core(TM) i7-7700 @ 3.6 GHz computer.

To underline the performance of the SNCR, a parametric study is then performed on the intentionally 435 mistuned structure presented in Section 5.2.2. The maximum displacements over all the sectors are represented in Figure 13 for differents values of excitation amplitude (ranging from 0.01 to 0.25 N). Such simulation lasted approximatively a day for the SNCR but would have required several weeks for the reference method.

In Figure 13, the maximal displacement over all sectors is represented as a function of the excitation frequency and the amplitude excitation for a linear and a nonlinear simulations. Using these results, different observations can be made to characterize the nonlinear system. It is for instance possible to determine the efficiency of the nonlinearity all across the excitation amplitude range. For the highest excitation amplitude simulated, the nonlinear effects assure more than a 50% displacement decrease at the resonant peaks.

Moreover, the nonlinearity only gets distinctly activated for an amplitude higher than 0.1 N. Notice also that standard features of the nonlinear system can be observed in Figure 13b: the shape of the backbone curve is similar to a quadratic shape (against a linear variation for the stuck case in Figure 13a), and the resonant peaks get flattened and widen by the nonlinearity activation level.

This result can now be compared to the nonlinear tuned case (system with only sectors of type 1) whose results are illustrated in Figure 14. To highlight the frequency split incoming from the mistuning effect, the results represent the displacement of the 8 th sector in the logarithm scale. As expected, a single ray is observed for the tuned case (see Figure 14a), whereas the mistuned system contains multiple rays (three in this case, see Figure 14b). Moreover, we can see, for this sector and the mode under study that the displacement amplitude is lower for the mistuned case than its tuned counterpart.

Through this short analysis, we can highlight that the SNCR is a suited candidate to tackle a wide range of applications for nonlinear cyclic structures and is able to perform accurate parametric studies which were, up to now, not realistically feasible on standard computers.

Conclusion

This paper proposes a new reduced-order model to determine the steady-state vibration of nonlinear mistuned structures. A synthesis procedure based on the cyclic complex nonlinear mode was first proposed for cyclic symmetric structures. This method was shown to be very accurate but only holds for traveling wave excitations. Combining a substructure technique applied to the system with the cyclic complex nonlinear mode, a new ROM named Substructuring method based on Nonlinear Cyclic Reduction (SNCR) was proposed. For all cases simulated, the new method was shown to be very accurate (below 2% of error) while reducing significantly the computation time (overall 15 times lower than the reference method). It was also compared with a recent ROM and was proved to be a more suitable candidate for nonlinear mistuned structures (error divided by 7 for a standard traveling wave excitation). The SNCR employs different assumptions, the main ones of which are the creation of a nonlinear basis composed of only specific cyclic NNMs and the substitution of the nonlinear forces by the formulation of the NNMs.

The flexibility of the SNCR is supposed to hold for a wide range of studies, applications and complex phenomena: underplatform dampers, gyroscopic effects, veering area, etc. Its reduced computation

  M . The degrees of freedom (DOFs) that are kept as master DOFs u master,1 are: the right boundary DOFs (subscript b,R), the left boundary DOFs (subscript b,L) contained in u b = u T b,R , u T b,L T , the nonlinear DOFs of both the blade and the corresponding portion of the disk (noted u nl ), and an observation DOF (u tip ) located at the tip of the blade. The internal modes of vibration are noted η and their number is chosen sufficiently high to ensure convergence of the associated linearized reduced model (considered as the stuck bladed-disk).

Figure 1 :

 1 Figure 1: Description of the reference and the CNCMS methodologies.

14090 × 7 =

 7 15120 unknowns (24 sectors, 90 relative nonlinear DOFs per sector and 3 harmonics).

  Fundamental sector with boundary DOFs. u nl,blade u nl,disk (b) Description of the contact interfaces.

Figure 2 :

 2 Figure 2: Description of the fundamental sector. The entire structure is made of 24 sectors.

Figure 3 :

 3 Figure 3: Frequency forced response ( ): linear response; ( ): nonlinear reference solution; ( ): CNCMS solution.

Figure 4 :

 4 Figure 4: Frequency forced response with a forward traveling excitation of different amplitudes. ( ): backbone curve of the cyclic complex nonlinear normal mode; ( ): nonlinear reference solution; ( ): cyclic synthesis method

Figure 5 : 4 . 1 .

 541 Figure 5: Frequency forced response for all sectors with a standing wave excitation (amplitude 0.25 N and h ex = 3). ( ): backbone curve of the cyclic complex nonlinear normal mode; ( ): nonlinear reference solution; ( ): cyclic synthesis method

Figure 6 :

 6 Figure 6: Schematic illustration of the SNCR.

Frequency

  (Hz) Amplitude (m) (b) Basis {0, 1, 3}.

Figure 7 :

 7 Figure 7: Frequency forced responses for all sectors under a standing wave excitation and h ex = 1. ( ): linear reference response; ( ): nonlinear reference solution; ( ): SNCR solution.

Figure 8 :

 8 Figure 8: Frequency forced responses for all sectors under a standing wave excitation and h ex = 3. ( ): linear reference response; ( ): nonlinear reference solution; ( ): SNCR solution.

Figure 9 :

 9 Figure 9: Frequency forced responses for all sectors under a localized excitation with an amplitude of 0.125 N. ( ): reference solution; ( ): SNCR solution.

  Linear response (stuck system).

  Zoom on Figure 10a (linear response).

  Zoom on Figure 10c (nonlinear response).

Figure 10 :

 10 Figure 10: Frequency forced responses for all sectors with random mistuning (σ = 15%).( ): full model solution; ( ): linear reference solution, ( ): SNCR solution, ( ): tuned response.

Figure 11 :

 11 Figure 11: Frequency forced responses for different values of σ . ( ): linear response for the reference solution; ( ): nonlinear reference solution; ( ): CNCMS solution; ( ): linear SNCR response; ( ): nonlinear SNCR solution.

Figure 11 representsFigure 12 :L

 1112 Figure 11 represents the maximum displacement of the entire structure as a function of the excitation frequency obtained with random mistuning defined with a standard deviation equal to 5%, 10% and 15%.Once again the SNCR provides a maximum amplitude very close to the reference solution for both the

Figure 2

 2 Figure 2 and 8b Figure 9b Figure 11c Figure 12ReferenceNo. of unknowns 15120 15120 15120 15120

Figure 13 :

 13 Figure 13: Parametric study on the amplitude excitation for a standing wave (h ex = 3).

Figure 14 :

 14 Figure 14: Comparison between the tuned and mistuned nonlinear displacement of the 8 th sector in the logarithm scale.

Table 1 .

 1

	State	Full	Reference	CNCMS
	States 0-1 40 internal modes 40 internal modes	40 internal modes
	States 2-3	-	600 internal modes 1 NNM with fixed boundaries
	States 3-4	-	-	400 interface modes

Table 1 :

 1 Components of the reduction basis for the full, the reference and the CNCMS methodologies.

Table 4 :

 4 Performance of the different methodologies on the simulations performed in Section 5. The simulations were run on a Intel(R)

The development presented in[START_REF] Quaegebeur | Model reduction of nonlinear cyclic structures based on their cyclic symmetric properties[END_REF] is based on a theoretical result which determines the nodal diameters coupled by friction nonlinearities and multiple wave numbers excitations.

Acknowledgements

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

linear and nonlinear cases. For better readability, these results only represent the maximum response of the entire structure as a function of the excitation frequency (and not each blade displacement separately). The main results of Figure 10 are thus summarized in Figure 11c in which the CNCMS results are also given for comparison. Similarly to what was shown in Figure 3, the CNCMS is not able to retrieve the correct maximal displacements.

Intentional mistuning

Some studies, such as the one in [START_REF] Castanier | Using intentional Mistuning in the design of turbomachinery rotors[END_REF], have shown that creating an intentional mistuning pattern may reduce the amplification factor. Intentional mistuning consists in assembling sectors which have been designed with different geometries or are made with different materials. In this study, two kinds of sectors are considered (M = 2) arranged with the following pattern taken randomly:

The number 1 denotes the standard sector used in Sections 2.5, 3.2 and 5.1, whereas 2 is a new sector derived from Sector 1 in which the length of the blade has been reduced by 0.5 mm. The reference system is obtained following the methodology given in Section 2 with m 1 = 14 and m 2 = 10. For the SNCR, the cyclic matrices of the two types of sectors are evaluated, and their respective NNMs are then computed. For these two kinds of sectors (i ∈ 1, 2 ), the basis Φ i in Equation ( 41) is formed by the NNMs corresponding to the sector of type i. In this test case, the system is excited by a standing wave excitation with h ex = 3 and an amplitude of 0.25 N. Hence, Φ i contains the NNMs of the zeroth, third and ninth nodal diameter of the sector type i (as explained in Section 4.2). Figure 12 represents the frequency forced response of the intentionally mistuned system. Once again, the results obtained with the SNCR are close to those obtained with the reference method. This validates the new methodology for intentionally mistuned structures.

and,

The quantity L ∞,u denotes the absolute relative error between the maximum displacements while L ∞,ω corresponds to the absolute relative error between their associated frequencies. ME u represents the average error over all the sectors and ME ω gives the average difference for their associated frequencies. Table 4 provides these information as well as the computation time and the number of unknowns for the nonlinear reference method, the SNCR and the CNCMS. The CNCMS performances are provided in Table 4 even though their associated results have not always been illustrated.

For all the test cases, the SNCR shows excellent attributes: its accuracy is very good and its computation time is at least 15 times slower than the reference solution. The highest error of 2% is reached when localization phenomena occurs (due to a localized external force). This can be explained by the fact that both mistuning and nonlinear phenomena couple many nodal diameters. As a consequence, even if the linear solution is converged, the nonlinearities may create additional couplings which will not be correctly captured by the reduction basis. However, these discrepancies are overall very small. The memory space required by the SNCR is also interesting as only the data for the M different sectors have to be saved during the solution procedure. In comparison, the CNCMS must store the model of the full disk with the boundary DOFs to perform the interface modes computation.

Overall the SNCR provides more accurate results and requires less computation time than the CNCMS method. As an example, for the simulation in Figure 8b the error level is reduced from 8.2% to 1.1% using the SNCR and the computation time is divided by 2.6. The only scenario where the CNCMS behaves better is for the localized excitation. In this particular case, the energy is mostly contained within a single sector and hence the basis of the CNCMS (fixed-boundary NNM + static modeshapes) is more accurate than the basis of the SNCR (composed of cyclic NNMs).

The main drawback of the SNCR is the computation time of the NNMs (around 50 min per NNM), as shown in Table 3. This is mainly due to the inversion of the linear dynamic rigidity matrix (containing the boundary DOFs) for the Schur condensation.

time makes it a very good candidate for nonlinear parametric studies to predict, for instance, the impact of mistuning and excitation amplitude on the amplification factor. The method is therefore expected to be of high interest for the scientific community working in the field of cyclic structures (perfectly symmetric or mistuned).