Arnaud Sternchüss
email: arnaud.sternchuss@ecp.fr

Keywords: structural dynamics, vibrations, turbomachinery, bladed disks, cyclic symmetry, mistuning, model reduction, substructuring, superelements

Bladed disks found in turbomachines are complex structures whose vibration characteristics are generally determined by exploiting the symmetry properties of their nominal configuration. This symmetry no longer exists either when disks are assembled to form a rotor or when discrepancies in the mechanical parameters are introduced (intentional or unintentional mistuning). Fine meshes required to correctly evaluate stress distributions would lead to prohibitive model sizes (typically a few million degrees of freedom). The objective of this thesis is to introduce model reduction techniques that rely on the combination of separate computations of acceptable size. This provides a means for in-depth studies of the behaviour of dense 3D models of multi-stage bladed rotors with possible mistuning. At first, Fourier transforms performed separately on each individual disk allows to understand the inter-harmonic coupling induced by inter-stage coupling and mistuning. From this study, a first method uses cyclically symmetric solutions plus sector modes with fixed inter-sector interfaces to build a reduced sector model. The latter is exact for target modes and very accurate for others. This method is extended to multi-stage assemblies by employing multi-stage mono-harmonic eigensolutions. Illustrations focus on the proposed methodology that enables to deal with large scale industrial models while remaining compatible with various post-processing procedures (free or forced response computations, analysis of their spatial harmonic content, energy distributions and localization effects...). This methodology is finally extended to the handling of parametric models depending on the rotation speed. The enrichment of the initial sets of target vectors with computations at three rotation speeds enables a fast and accurate recovery of the evolution of the eigenfrequencies with respect to the rotation speed in any operating range.

Résumé

Les disques aubagés, que l'on trouve dans les turbomachines, sont des structures complexes dont le comportement vibratoire est généralement déterminé par l'exploitation de conditions de symétrie dans leur configuration nominale. Cette symétrie disparaît lorsque l'on assemble plusieurs de ces disques pour former un rotor ou que l'on introduit une variabilité spatiale des paramètres mécaniques (on parle de désaccordage intentionnel ou non). Le raffinement des maillages, nécessaire à une évaluation correcte de la répartition des contraintes, conduirait à des modèles de rotor complet de taille prohibitive (plusieurs dizaines de millions de degrés de liberté). L'objectif de cette thèse est donc l'introduction de méthodologies de réduction qui par combinaison de calculs acceptables permettent d'étudier de façon fine la dynamique d'ensemble sur des modèles 3D fins multi-étages et potentiellement désaccordés. L'étude des transformations de Fourier séparées des réponses de chaque étage permet, dans un premier temps, de bien comprendre les effets de couplage inter-harmonique liés au couplage inter-disque et au désaccordage. A partir de ce constat, une première méthode utilise les résultats de calculs en symétrie cyclique et à secteur encastré pour construire un modèle de secteur exact pour certains modes dits cibles et de très bonne qualité pour les autres modes. Cette méthode est ensuite étendue au cas multi-étage en construisant des bases de réduction de secteur par combinaison de solutions mono-harmoniques. Les illustrations montrent que la méthodologie proposée permet le traitement de modèles de très grande taille, tout en restant compatible avec une grande richesse de post-traitements (calculs de modes, calculs de réponses forcées, analyses de leur contenu harmonique spatial, répartition d'énergie et effets de localisation...). La méthodologie est enfin étendue à la gestion de modèles paramétrés en vitesse de rotation. L'enrichissement des ensembles de modes cibles par des calculs à trois vitesses permet ainsi une reconstruction rapide de l'évolution des fréquences pour l'ensemble d'un intervalle.

Mots-clés : dynamique des structures, vibrations, turbomachines, disques aubagés, symétrie cyclique, désaccordage, réduction de modèle, sous-structuration, super-éléments.

i

Introduction General description of turbomachines

A gas turbine is a manufactured machine that exploits the mechanical work created by the continuous combustion of an air/fuel mix [Rol 2005]. Air is swallowed at the intake of the engine by the fan. In engines without bypass, all the air enters the core of the engine (Figs. 1c to e). Otherwise, a fraction of the air flow bypasses the core of the engine (Fig. 1a andb). The core flow passes through several stages of the low pressure compressor (LPC). On three-shaft engines, the fan is the only stage of the LPC and the air subsequently goes through an intermediate pressure compressor (IPC). This flow is compressed further in the high pressure compressor (HPC) before it is mixed with fuel. This mix burns in the combustion chamber and hot gases are then expanded in the high pressure turbine (HPT), the intermediate pressure turbine (IPT), if any, and the low pressure turbine (LPT). Each compressor is driven by its own turbine stage(s) which gives a control of its angular speed. From high to low angular speeds one finds the HPC/HPT, the IPC/IPT, if any, and the LPC/LPT.

The power transmitted to the driving shafts is only a small fraction of the power produced by the engine. The rest of this power is used for various applications:

• in aviation, it takes the form of a thrust that pushes the engine and the attached plane forward in turbofans with either a high bypass (Fig. 1a) or low bypass (Fig. 1b) ratio and turbojets, a torque that drives a shaft attached to a propeller that produces thrust in turboprops (Fig. 1c), a torque that drives both the main rotor, which produces lift and thrust, and the auxiliary rotor, that ensures stability in helicopter turbines (Fig. 1d).

• in marine applications, it takes the form of a torque that drives either a propeller or a water turbine, • in ground applications, it takes the form of a torque that drives a generator in power plants (Fig. 1e), a gas/oil pump in pipeline applications or a power shaft in military vehicles.

For all applications, a turbomachine is made of a static casing attached to the vessel or the ground inside which two or three shafts are rotating concentrically (see Fig. 2a), with contra-rotation in some engines [Rol 2005]. In the engine displayed in Fig. 2a, the inner shaft connects the fan and the LPC to the LPT. The outer shaft connects the HPC to the HPT. In compressors and turbines, the air flows through interlaced rotor and cantilevered or mobile stator vanes.

A compressor (resp. turbine) stage is made of a pair (rotor disk,stator disk). In a compressor, the rotating blades impart energy to the air resulting in an increase in its static pressure. However, this operation swirls the flow and the latter has to be deswirled by stator vanes so that it enters straight into the following compressor or turbine stage. Variable inlet guide vanes and variable stator guide vanes are used to control the airflow in the compressor. The pressure ratio of a single stage is quite small, 1.2:1 ¡ 1.3:1, but the high number of stages enables to reach an overall ratio of 30:1 ¡ 40:1 (see Fig. 2b). The inverse phenomenon occurs in a turbine stage since the power is extracted from the air leading to a decrease in the static pressure. Note finally that the bypass flow in a turbofan may contribute up to 70% of the net thrust [Rol 2005]. A rotor is made of a compressor drum connected to a turbine disk or a turbine drum. A compressor drum is an assembly of bladed disks with commonly 3 to 12 stages. These disks are bolted or welded together so that it forms a fairly stiff assembly. Seal annuli are machined in the inter-stage rims so that a sufficient tightness is ensured at the tip of the matching stator blades.

A compressor or turbine bladed disk is made of an inner cylindrical part (the disk) and an outer bladed part. Different technologies exist:

• blades can be mounted on a disk machined from a separate cylinder of metal. Two configurations exist: blades can be either axially mounted (first stages) or circumferentially mounted (last stages). In both cases, securing devices ensure that the blades remain safely mounted during operation. The advantage of such a configuration is that blades can be individually chosen for a given reason (limit in the discrepancies due to machining tolerances between adjacent blades) or individually replaced in case of failure. This also allows the design of wide-chord hollow composite blades for the fan which are lighter than their metallic counterpart but at least as resistant. In turbines, blades are micro-perforated so that air supply can be used to cool them down, otherwise they would operate beyond their melting point even with a thermal barrier coating. Finally, the friction interfaces provide significant damping to the disks for energy dissipation. However, this solution is relatively heavy due to the presence of the securing devices.

• integrally bladed disks (blisks) are machined from a single block of metal (Fig. 3). The enormous advantage of this solution is that it provides lighter disks. It is however very sensitive to the machining tolerances and discrepancies in the geometry from one blade to another are unavoidable, especially because of tool wear during the machining process [Feiner and Griffin 2004b]. Besides, these disks are very lightly damped which is a key issue in vibration.

Figure 3: Blisk in a GE90-115B

• integrally bladed rings (blings, currently in development) further reduce the overall weight since the inner rim is replaced by a stiff ring. Blings would also be lightly damped by construction.

Jet and helicopter engines have very specific design requirements. Since their overall mass is a key parameter in the fuel consumption of the flying vehicle they are mounted on, their components have to be as light as possible. This is why light material with high resistance such as titanium alloys are employed for the cold parts. As the temperature increases in the engine (see Fig. 2b), steel alloys coated with heat-resistant materials, such as ceramics, are used especially in the turbine stages. Moreover, some parts are very slender which may make them very soft and sensitive to the inertial loading.

Design constraints

The prediction of mechanical responses of bladed rotors requires a good knowledge of the external loading to which they are subject. Details can be found in [Rol 2005, Srinivasan 1997, Slater et al. 1999, Castanier and Pierre 2006]. The rotating parts have to withstand:

• static loads:

inertial loads induced by the rotation of a heavy structure with a high angular speed which tends to stiffen the structure as the rotation speed increases,

loads due to a global thermal gradient between the turbine and the front stages of the compressor (see Fig. 2b) and a local gradient from the blade to the inner rim of the disk which tends to soften the structure as the temperature increases, weight and steady pressure loads that are typically negligible compared to thermal and inertial loads;

• to which dynamic loads are superimposed:

high aerodynamic forces (stationary or non-stationary),

bearing loads and inter-stage loads, contact loads when a blade tip impacts the casing, impact loads if a foreign object goes through the air flow.

The aerodynamic forces are of specific interest since they are induced by the upstream and downstream flow due to the presence of equally spaced obstacles such as the stator or nozzle guide vanes. As a result, the aerodynamic loading is classically computed with a computational fluid dynamics (CFD) code in cyclic symmetry. The amplitude of the excitation applied to a rotor disk is distributed circumferentially according to a cosine pattern whose angular frequency δ is the number of obstacles right upstream to this disk. In the frame attached to the rotating disk, it is a so-called backward-travelling engine order excitation whose frequency is an integer multiple of the rotation speed since a blade sees δ excitations per revolution [START_REF] Lalanne | [END_REF]. One reports in a Campbell diagram the evolution of the eigenfrequencies with respect to the angular speed Ω and the engine order lines, i.e. the straight lines ω δΩ. The abscissa where two curves in these sets intersect defines a possible resonant speed. Otherwise, the distance between the eigenfrequency curves and the engine order lines defines a frequency safety margin. One contribution of the present study is to provide an accurate technique to build the Campbell diagrams so that resonances can be tracked.

Therefore, one of the main causes of turbomachine failure is the individual blade failure due to High-Cycle Fatigue (HCF). When a resonance occurs on a Campbell diagram, blades withstand a certain number of cycles of alternate stresses (see small graph in Fig. 5). The length of time needed to reach a high number of cycles depends on the resonance frequency, it can be only a few seconds if this frequency is high enough. Fatigue theory predicts the maximum level of dynamic stress σ D the blades can withstand during N c cycles when the static (or mean) stress σ S is known. When the stress state can be considered uniaxial, the Haigh diagram gives the safe region in which the integrity of the blades is certain, e.g. no less than 99.99% in a blade population will survive. The boundary of the safe region is approximated by a certain number of curves displayed in Fig. 5. A pair (σ S , σ D ) is attached to each operating point. For a given σ S , the distance of σ D to the chosen boundary curve defines the stress safety margin. If the margin

Ω Ω max ω j 1Ω 2Ω 3Ω 4Ω 5Ω
Figure 4: Campbell diagram with (--) engine-order excitation lines is zero or negative the blade geometry has to be updated so that it satisfies this constraint. The key point in HCF analysis is to compute accurately both σ S and σ D . Multi-axial fatigue is somewhat more complex since it involves the Von Mises or Tresca equivalent stresses in much more complex criteria. Damping also plays a key role in HCF. For low damping the amplitudes, hence the stresses, are inversely proportional to damping. Thus special attention must be paid to this parameter for integrally bladed disks since their structural damping is very low, in contrast with disks with inserted blades that present a "natural" dissipation by friction at the contact interfaces [Petrov and Ewins 2006]. In both cases, the overall damping is artificially increased by friction devices commonly mounted at the root of the blades, under the platforms [Petrov 2007, Sever et al. 2007], in grooves machined in the inter-stage rims of blisks [Laxalde et al. 2007b, Laxalde et al. 2007c,Laxalde et al. 2008] or added as snubbers and shrouds located between half the length (compressors) and the tip (turbines) of the aerofoils [Chen andMenq 2001,Petrov and[START_REF] Petrov | [END_REF]. Visco-elastic materials and piezo-electric materials are also good candidates for future applications.

Other causes of failure are beyond the scope of this work. From [Rol 2005] one can cite for the sake of completeness:

• rotating stall which occurs at low speed. Stall is a sudden decrease of the lift of the blade when its critical angle of attack is exceeded. This phenomenon may arise due to an excessive amplitude of vibration of the aerofoil which acts on the angle of attack. Part span stall considers a few unconnected stall cells, full span stall considers a single cell spreading over a wide area of the rotor. Rotating stall affects mostly the first stages in a compressor, it gives rise to a flow that is circumferentially non-uniform and rotates at 20 to 50% of the angular speed of the rotor.

• surge which occurs at higher speeds. Surge refers to the oscillation of the air flow in the compressor between stalled and unstalled flow. The loss of pressure between the intake and the outtake of the HP compressor may even lead to the forward ejection of the combustion gases.

• flutter which also occurs at higher speeds, when the aerodynamic flow becomes unsteady.

Flutter is a self-excited oscillation at a frequency close to that of the blades. As a result, blades may encounter increasing levels of stress resulting in their rapid failure [Srinivasan 1997]. Flutter may be seen as a negative damping since energy is continuously provided to the blades by the means of the fluid-structure coupling terms. The instability is set when the structural damping cannot compensate for this negative aerodynamic damping.

• foreign object damage which occurs when either the low compressor blades are impacted by a bird or a block of ice at takeoff or landing, or if a downstream blade is impacted by debris of an upstream blade. Whatever the case, the metal pieces must be contained by the casing. Fan blade-off and bird ingestion are two tests used in the certification process of an engine.

As a conclusion, the different tasks required for failure analysis in an industrial design process can be summarized in the flowchart of Fig. 6.

Motivations

The present study is only interested in the left-hand path of this flowchart. One thus needs to accurately compute resonance frequencies for Campbell diagram margins and stress distributions for fatigue predictions.

Current procedures use cyclic symmetry to predict responses of individual disks/stages. To account for inter-disk coupling, axisymmetric models of the full rotor give access to static deformations that are fairly accurate in the rims and can be used as enforced motion of a single disk/stage. When computing modes, a single disk is however considered so that no inter-disk coupling can be predicted.

In practice, newer designs with more flexible rims lead to cases where not all modes are strongly confined to a single stage. One can thus have high responses on a stage due to excitation on another. The ability to predict such coupling and the analysis of its occurrences on realistic models is the first motivation for the present work. Such computations cannot be performed using standard cyclic symmetry since the number of blades in each stage is typically different. A first contribution of the present work is therefore to present a consistent way to perform Fourier decompositions of the responses independently for each stage. The idea of computing spatially mono-harmonic responses [Laxalde et al. 2007a] then appears as a simple extension of Figure 6: Flowchart of the structural design process cyclic symmetry to multi-stage applications. There are however underlying assumptions that are detailed further.

When considering meshes that are sufficiently refined for fatigue analysis, the sector mesh typically uses 10 5 ¡ 10 6 degrees of freedom (DOF). Without cyclic symmetry, a full rotor model would use many millions of DOF for highly connected volume models. The perspective to practically run such models is very remote. A second contribution is consequently the introduction of a reduction methodology leading to accurate predictions of the full order 3D model of the full rotor by combining the results of a series of problems that can currently be solved.

In classical Component Mode Synthesis (CMS) methods, components are assumed to be reduced independently. This leads to a reduction as a series of decoupled problems but also induces a loss of accuracy. A first change proposed here is to use target mono-harmonic solutions as initial Ritz vectors for the reduction.

In mono-disk analyses, this enables exact tuned responses for retained modes and very accurate predictions of all modes for both tuned and mistuned applications. The ability to select target modes is seen as a very useful advance for industrial applications, where the target loads are typically known.

A second change is to consider disjoint models connected by an interface without internal DOF. This change of perspective significantly facilitates the expression of coupled responses. For multi-stage rotors, mono-harmonic rotor modes are thus combined to lead to reduced sector models that are independent for each stage. The methodology is shown to be practical for both academic test cases and an industrial rotor model, with second order tetra10 elements and a total of 22 million nodes.

Beyond the results needed to demonstrate the accuracy of the proposed methods, the focus is set on the needs for a detailed methodology to apply such methods to large scale industrial models. In particular, specific post-processing procedures have been set up to handle recoveries of the modeshapes to various post-processing meshes.

The last contribution is linked to the ability to rapidly predict frequencies at multiple rotation speeds. In applications with slender parts, the geometrical stiffening linked to pre-stress can induce significant frequency shifts and modeshape changes as speed varies. A multi-model reduction approach is thus introduced to allow very accurate predictions of modes and frequencies as a function of rotation speed based on computations at only a few speeds.

Manuscript Outline

The nominal geometry of a bladed disk presents a very convenient property of geometrical periodicity called cyclic symmetry, where the full geometry derives from that of a single angular sector. This condition is sufficient to treat the mechanical problems in the framework of the spatial Discrete Fourier Transform (DFT). Chapter 1 presents the theoretical results associated with such a treatment and focuses on the convenience of the DFT when seen as a matrix transformation in all the proposed applications. When further restricting the problem to the tuned case where the mechanical properties (mass, stiffness, damping,...) are identical from a sector to another, one demonstrates the classical result that the full matrix problem, once projected, can be subdivided into independent sub-problems attached to the involved Fourier harmonics of the transform of displacements and external forces. A classical consequence of this transformation is that tuned bladed disks have spatially mono-harmonic modeshapes. When the properties are not periodic, the disk is said to be mistuned. The Fourier harmonics are shown to be coupled leading to modeshapes with a rich spatial harmonic content. Rotors are then introduced as assemblies of cyclically symmetric disks. Except in very particular cases, the resulting structure does not belong to any group of symmetry so that the finite element matrices of the rotor do not present any symmetry. Considering that the geometry of each disk is periodic, Fourier transforms of the associated DOF can be performed leading to decoupled problems for each disk coupled by a fully populated interface matrix that describes the energy of the rim connecting the disks. Ignoring inter-harmonic coupling at the interface, leads to families of mono-harmonic rotor problems which can be solved for eigenvalues. The validity of the assumption that inter-harmonic coupling can be neglected is a key issue addressed in the following chapter.

Chapter 2 enters into the details of the two proposed reduction techniques suited to the study of single-and multi-stage bladed rotors. The two methods rely on the same type of substructuring: each disk is subdivided into sectors to which superelements are attached. The kinematic subspace in which the motion of the superelement is sought is spanned by a set of vectors that arises from the transformation of spatially mono-harmonic eigenmodes (resp. eigensolutions) computed with the single stage (resp. multi-stage) cyclic symmetry approach plus sector modes with inter-sector (resp. inter-sector/inter-disk) interfaces fixed. The resulting reduced model converges towards the true motion of the structure with an accuracy controlled by the initial set of target vectors. The multi-stage methodology is fully developed allowing predictions of modes and forced responses, full and partial recovery of solutions. This is used to discuss issues in the possible applications. Restitution to display meshes enables a straightforward check of the deformed shapes or general forced response computations. Restitution to partial meshes is used for stress analysis. Direct processing of the reduced model generalized DOF is used to analyze harmonic content and disk contributions.

Finally, Chap. 3 presents the multi-model reduction technique to compute the response of a rotating structure for a range of rotation speeds. The first part of this chapter paves the way for the approximation of the stiffness matrix by a matrix polynomial by recalling how the elastic stiffness is supplemented with centrifugal and geometric terms as the rotation speed increases. It is then shown that three sets of exact solutions computed at target rotation speeds are sufficient to accurately represent modes and frequencies at all intermediate values. A typical application is presented where a subspace spanned by the target vectors is used as the input of a subsequent reduction process. It leads to well-conditioned problems with the help from the same orthogonalization procedures developed in the previous chapter. The accuracy of the multi-model approach is demonstrated for two different reductions techniques. These examples provide the occasion to describe the associated methodology. The present chapter proposes a unified presentation of cyclic symmetry for single disks and multi-stage assemblies. The base assumption is that the geometry is rotationally periodic. In other words, each stage can be divided into angular sectors that are generated by successive rotations of the initial (or generating) sector. Under this assumption, one can perform spatial Discrete Fourier Transforms on the equations, which lead to simplifications detailed in the chapter. Since the geometry is periodic, one can define series of displacements at nodes corresponding to the periodicity. Rather than using the series of physical DOF, one can use the coefficients of its Discrete Fourier transform [Ewins 1973, Henry and Ferraris 1984, Ohayon and Soize 1998]. This transform is described in certain papers as a Floquet expansion [Touratier 1986, Lalanne and[START_REF] Lalanne | [END_REF]] or a result of the finite group theory [START_REF] Capiez-Lernout | Dynamique des structures tournantes à symétrie cyclique en présence d'incertitudes aléatoires[END_REF], but it is basically the same theory.

1

Dynamics of Geometrically Periodic Structures

Sectors are first assumed to share the same set of mechanical properties so that the resulting structure is said to be tuned. In that case, thanks to the uncoupling of the Fourier harmonics, the model equations are decoupled by block and can be related to a matrix problem that involves a one-or two-sector approach with the proper inter-sector dephasing conditions [Thomas 1974, Thomas 1979,Henry and Ferraris 1984,Jacquet-Richardet et al. 1996,Yang and Griffin 1997[START_REF] Capiez-Lernout | Dynamique des structures tournantes à symétrie cyclique en présence d'incertitudes aléatoires[END_REF], Chatelet et al. 2005]. Such a decomposition is commonly found in industrial codes [START_REF] Msc | MSC Software[END_REF], Nicolas et al. 2006[START_REF] Sam ; Samtech | Rotor module for the analysis of rotating systems[END_REF].

Section 1.1 recalls this classical theory with a particular matrix formulation suited to the later study of multi-stage assemblies of bladed disks. Extensions not considered here are methods that involve Lagrange multipliers [Irwanto et al. 2001] or methods with non-linear effects (shrouds, snubbers, friction dampers...) [Chen and Menq 2001, Petrov and Ewins 2005, Nikolic et al. 2007, Petrov 2007, Sever et al. 2007, Laxalde et al. 2007b, Laxalde et al. 2007c, Petrov 2008].

In Sec. 1.2, one is interested in mistuned structures, that is to say structures whose mechanical properties vary from a sector to another. In this case, it is no longer possible to decouple the Fourier harmonics and full finite element models are required, except for very specific mistuning configurations described in that section. Mistuning corresponding to changes in geometry is not considered since the assumption of geometric periodicity is not verified.

Section 1.3 addresses multi-stage assemblies of tuned bladed disks. The main part of this section is devoted to disks connected with a volumic ring that concentrates all the asymmetry of the system. The concept of Multi-Stage Cyclic Symmetry (MSCS) introduced in [Laxalde et al. 2006, Laxalde et al. 2007a] is recalled in that specific framework. It allows to define mono-harmonic eigensolutions (by analogy to the mono-harmonic modes of a single tuned disk) computed from a matrix problem that is similar to that used for a single structure. Such eigensolutions, if they are not eigenmodes, are however good approximations of them.

Fourier Decomposition on a single periodic geometry

Fourier transform of fields

A structure is said to present cyclic symmetry when it is composed of N s geometrically identical angular sectors, numbered from 0 to N s ¡ 1, that are generated by rotations of a reference sector S 0 whose angular aperture is α 2 π{N s such as that of Fig. 1.1. One insists that the only requirement of this section is that sectors are geometrically identical even though they may have different mechanical properties. x 0 stands for the position of any point in sector S 0 . Let θ stand for the rotation of angle α around axis Oe z . The application of successive rotations θ s for s in v0, N s ¡ 1w to x 0 then generates a series px s q, s in v0, N s ¡ 1w, such that x s belongs to sector S s . The values of any field v defined on the full set of points px s q form a finite series pvpx s qq. As described in Appendix B.1, one can compute the unitary discrete Fourier transform of each component of this field, thus leading to the unitary discrete Fourier transform of the full field which is given by

p v 0 px 0 q 1 c N s Ns¡1 ş0 vpx s q, p v Ns 2 px 0 q 1 c N s Ns¡1 ş0 vpx s q p¡1q s , (1.1)
and other harmonics are such that dδ v1,

N s {2 ¡ 1w, p v δ px 0 q 2 c N s Ns¡1 ş0
vpx s q pcosps δ αq i sinps δ αqq .

(1.

2)

The fields p v δ px 0 q are called the Fourier harmonics, they are associated with a Fourier harmonic coefficient δ. The series pvpx s qq is recovered by the unitary inverse Fourier transform, ds v0, N s ¡ 1w,

vpx s q 1 c N s p v 0 px 0 q Ns 2 ¡1 δ1 2 N s Re p v δ px 0 q ¨cosps δ αq ¡ Im p v δ px 0 q ¨sinps δ αq ¨ p¡1q s c N s p v Ns 2 px 0 q. (1.3)
For any s in v0, N s ¡1w, vpx s q is a sum of the Fourier harmonics p v δ px 0 q multiplied by a dephasing factor that depends on the value of the Fourier coefficient δ and the sector number s. In other words, two adjacent sectors vibrate with an angular phase difference of δ α [Srinivasan 1997,Pierre and[START_REF] Pierre | [END_REF].

Notice that this transform can be applied to any field defined on the geometrically repetitive structure, including pressure, load, displacement, stresses, strains... In the following, the matrix representation introduced in Appendix B. 1.4 for scalar functions will be preferably used. It concerns the coordinates of v taken separately. If one defines a vector tvu that contains all the terms in the series pvpx s qq and a vector tp vu that contains all the Fourier harmonics, that is to say if they are such that

tvu t 2 vpx 0 q ¤ ¤ ¤ vpx s q ¤ ¤ ¤ vpx Ns¡1 q @ , tp vu t 3 p v 0 px 0 q ¤ ¤ ¤ Re p v δ px 0 q ¨Im p v δ px 0 q ¨¤ ¤ ¤ p v Ns 2 px 0 q A . (1.4)
The unitary discrete Fourier transform of tvu is then given by tp vu

E ¡1 $ tvu , (1.5)
and the unitary inverse discrete Fourier transform of tp vu is tvu rEs tp vu .

( 1.6) rEs is the unitary Fourier matrix defined in Appendix B.1.4 whose expression is recalled here

rEs 1 c Ns " " " " " " ! 1 ¤ ¤ ¤ c 2 0 ¤ ¤ ¤ 1 . . . . . . . . . . . . 1 ¤ ¤ ¤ c 2 cos ps δ αq ¡ c 2 sin ps δ αq ¤ ¤ ¤ p¡1q s . . . . . . . . . . . . 1 ¤ ¤ ¤ c 2 cos ppN s ¡ 1q δ αq ¡ c 2 sin ppN s ¡ 1q δ αq ¤ ¤ ¤ p¡1q Ns¡1 ( 0 0 0 0 0 0 )
. (1.7)

In the rest of this chapter, rE s s stands for the block-row in rEs that corresponds to sector S s and rE δ s stands for the block-column in rEs that corresponds to Fourier harmonic coefficient δ.

Fourier transform of discretized fields

In finite element formulations (Ritz-Galerkin procedure), the continuous field v is discretized which leads to a physical vector of degrees of freedom (DOF) on the full structure. Since it is assumed that the geometry is rotationally periodic, one can group DOF as a series of vectors tv s u defined on a sector S s . Therefore, it becomes possible to compute the Fourier transform of the series of values at each DOF image of its counterpart in the generating sector. Let tvu be such that tvu t t v 0 v 1 ¤ ¤ ¤ v Ns¡1 u. It has the same discrete Fourier transform as the continuous field v which is

tp vu r p E I N 0 s tvu , (1.8)
where N 0 is the number of DOF of sector S 0 and tp vu is the harmonic DOF vector that can be defined such that

tp vu t 3 p v 0 Re pp v 1 q Im pp v 1 q ¤ ¤ ¤ Re ¡ p v Ns 2 ¡1 © Im ¡ p v Ns 2 ¡1 © p v Ns 2 A .
tvu can be recovered with the help from the inverse discrete Fourier transform tvu rE I N 0 s tp vu.

(1.9)

Notice that, at this point, the physical discretized field defined on two adjacent sectors tv s u and tv s 1 u are assumed disjoint. The continuity conditions are discussed in Sec. 1.1.4.

Fourier transform of matrices

The Fourier transform of a field, either continuous or discretized, only requires a periodic distribution of points. This section aims to examine the transformation of the matrix representation of the operators that act on these fields.

Let rA s s be any finite element matrix associated with sector S s . Matrix rAs results from the assembly of the sector matrices (recall that sector are still assumed disjoint) .10) In this section, one assumes that all sectors share exactly the same set of mechanical properties:

rAs " " " " " " ! A 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ A s ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ A Ns¡1 ( 0 0 0 0 0 0 ) . ( 1
the bladed disk is said to be tuned. In that particular case, matrix rAs derives from matrix rA 0 s associated with sector S 0 so that rAs .11) The congruent transformation of matrix rAs is then

I Ns A 0 $ . ( 1
r p As E t I N 0 $ I Ns A 0 $ rE I N 0 s E t E $ A 0 $ " " " " " ! A 0 0 ¤ ¤ ¤ 0 0 0 A 0 ¤ ¤ ¤ 0 0 . . . . . . . . . . . . . . . 0 0 ¤ ¤ ¤ A 0 0 0 0 ¤ ¤ ¤ 0 A 0 ( 0 0 0 0 0 )
.

(1.12)

The blocks associated with harmonic coefficients δ 0 and δ N s {2 are simply rA 0 s. The blocks associated with the other harmonic coefficients are

A 0 0 0 A 0 & ,
where the rows correspond to the real and the imaginary parts of the corresponding harmonic respectively.

To illustrate this, Fig. 1.2 describes a very simple lumped model that represents a single tuned disk. The initial sector is composed of several degrees of freedom. Its mass is concentrated at the first DOF and its first spring is two times stiffer than the others. If one considers a disk with 12 sectors and 4 springs per sector, the mass and stiffness matrices of its initial sector are M 0 $ " " " " ! m 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

( 0 0 0 0 ) , K 0 $ " " " " ! 2 k d ¡2 k d 0 0 0 ¡2 k d 3 k d ¡k d 0 0 0 ¡k d 2 k d ¡k d 0 0 0 ¡k d 2 k d ¡k d 0 0 0 ¡k d k d ( 0 0 0 0 )
.

(1.13) Figure 1.3 shows the projected mass and stiffness matrices

E t I N 0 $ rI N s M 0 s rE I N 0 s and E t I N 0 $ rI N s K 0 s rE I N 0 s.
The block-rows and block-columns correspond to the Fourier harmonics. For the sake of readability and further analyses, the color associated with a given harmonic δ starts at the block indexed by δδ and spreads towards the rows and the columns whose harmonic δ I is such that δ I δ. These projected matrices have exactly the shape described by Eqn. (1.12). The property exhibited in this section is very important since it states that if the parameters that drive the mechanical behaviour are identical from a sector to another, the matrices do not couple the Fourier harmonics to each other. It is equivalent to saying that subsequent harmonic coupling can be observed due to changes in the skeleton of the finite element matrices.

Continuity conditions between sectors

Let us consider that each sector S s can be decomposed into a complementary domain S s c and two matching right and left surfaces I s r and I s l so that for any x s r on I s r , x s l θpx s r q is on I s l , as displayed in Fig. 1.4.

An arbitrary field v defined on the full structure satisfies the inter-sector continuity conditions between sectors S s and S s 1 that are In the case when the field v is the displacement (resp. the load), a set of linear constraints associated with the interface DOF [Balmès 1996c] are given by the right and left observation equations which for compatible meshes give a discrete representation of upx s r q (resp. f px s r q) on I s r and θ pupx s r qq (resp. θ pfpx s r qq) on I s l such that upx s r q ùñ tq s r u rc r s tq s u , θ pupx s r qq ùñ tq s l u rc l s tq s u .

x s l x s 1 r θ px s r q ùñ vpx s l q vpx s 1 r q θ pvpx s r qq . ( 1
(1.15)

The rotation θ is included in rc l s so that the displacement is written in a basis coherent with the motion at the interfaces. If they are regularly meshed, the observation matrices are the same for all sectors. Thus, the discrete form of Eqn. (1.14) is given by rc l s tq s u ¡ rc r s 2 q s 1 @ t0u .

(1.16)

The constraint equation between right and left interfaces can be obtained in two ways:

(i.) since the Fourier harmonics are given by the discrete form of Eqns. (1.1) and (1.2), the restriction of the harmonics to each inter-sector interface can be computed which leads to dδ v0, N s {2w,

rc l s tp q δ u 2 c N s Ns¡1 ş0 rc l s tq s u e ¡i s δ α 2 c N s Ns¡1 ş0 rc r s 2 q s 1 @ e ¡i s δ α £ 2 c N s Ns¡1 ş0 rc r s 2 q s 1 @
e ¡i ps 1q δ α e i δ α rc r s tp q δ u e i δ α , which finally leads to dδ v0, N s {2w, rc l ¡ c r e i δ α s tp q δ u t0u , (1.17) or equivalently to dδ v0, N s {2w, rc r ¡ c l e ¡i δ α s tp q δ u t0u .

(1.18)

Once more, the harmonics are found to be uncoupled and a single constraint equation can be written separately for each harmonic such that rp c 0 s tp q 0 u t0u , rp c δ s 4 Re pp q δ q Im pp q δ q B t0u , rp c Ns 2 s tp q Ns 2 u t0u .

(1.19)

If one gathers the harmonics into the harmonic DOF vector, the global constraint equation is expressed as

rp cs tp qu t0u , (1.20)
which is exactly

" " " " " " " " ! c l ¡ c r ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ c l ¡ cospδ αq c r sinpδ αq c r ¤ ¤ ¤ 0 0 ¤ ¤ ¤ ¡ sinpδ αq c r c l ¡ cospδ αq c r ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ ¡pc l c r q ( 0 0 0 0 0 0 0 0 )
6 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 7 p q 0 . . . Re pp q δ q Im pp q δ q . . .

p q Ns 2 D G G G G G G G G F G G G G G G G G E
6 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 7 0 . . . .

0 0 . . . 0 D G G G G G G G G F G G G G G G G G E . ( 1
(1.23)

The substitution of the physical DOF vector tqu by its Fourier decomposition given in Eqn. (1.9) is rcs rE I N 0 s tp qu t0u .

(1.24)

It requires the computation of the product rcs rE I N 0 s, given by rcs rE

I N 0 s 1 c Ns rc l ¡ c r s ¤ ¤ ¤ 2 Ns rcosps δ αq ¡ sinps δ αqs c l ¡ c r cospδ αq c r sinpδ αq ¡c r sinpδ αq c l ¡ c r cospδ αq & ¤ ¤ ¤ p¡1q s c
Ns rc l c r s % rI N 0 s.

(1.25)

Equation (1.24) describes a system of linear equations where each line s, numbered from )), the lines can be properly summed to make the coupling between harmonics disappear:

0 to N s ¡
• the sum of all the lines multiplied by 1 leads to N s rc l ¡ c r s tp q 0 u t0u .

(1.26)

• the sum of all the lines multiplied by cosps δ αq ¡ sinps δ αq leads to dδ v1,

N s {2 ¡ 1w, 2 N s c l ¡ c r cospδ αq c r sinpδ αq ¡c r sinpδ αq c l ¡ c r cospδ αq & tp q δ u t0u .
(1.27)

A equivalent form of the latter can be exhibited since

c l ¡ c r cospδ αq c r sinpδ αq ¡c r sinpδ αq c l ¡ c r cospδ αq & cospδ αq sinpδ αq ¡ sinpδ αq cospδ αq & c r ¡ c l cospδ αq c l sinpδ αq ¡c l sinpδ αq c r ¡ c l cospδ αq & .
• the sum of all the lines multiplied by p¡1q s leads to

N s rc l c r s 3 p q Ns 2 A t0u .
(1.28)

As described in this section, the inter-sector continuity conditions do not combine harmonics together. They have to be added to the mechanical problem so that the vibration problem can be formulated for a single geometrically periodic structure.

Response of a tuned disk to external excitations

Let tqu be the discrete displacement, tfu the discrete external forcing and rZs the dynamic stiffness of a tuned bladed disk D with N s blades. In the undamped case, rZs rK ¡ ω 2 M s.

The problem on the full disk is rZs tqu tfu . (1.29) In this section, all sectors are assumed to have the same mechanical properties, thus rZs is obtained by rI Ns Z 0 s with rZ 0 s being the dynamic stiffness of sector S 0 (see Sec. 1.1.3) Two harmonic DOF vectors tp qu and t p f u are associated with tqu and tfu respectively when their Fourier transform is performed as described by Eqn. (1.8). Equation (1.29) is then turned into r p Zs tp qu t p f u, (1.30) which is exactly

" " " " " " " " ! Z 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ Z 0 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 Z 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ Z 0 ( 0 0 0 0 0 0 0 0 )
6 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 7 p q 0 . . . Re pp q δ q Im pp q δ q . . .

p q Ns 2 D G G G G G G G G F G G G G G G G G E 6
9 9 9 9 9 9 9 9 9 8 9 9 9 9 9 9 9 9 9 7 p f 0 . . .

Rep p f δ q Imp p f δ q . . . p f Ns 2 D G G G G G G G G G F G G G G G G G G G E . (1.31)
Inter-sector continuity is added to Eqn. (1.31) so that the full mechanical problem of Eqn. (1.30) can be decomposed into a class of parametrized problems that allow to compute the individual responses of the Fourier harmonics such that:

• if δ 0, tp q 0 u is real and satisfies Z 0 $ tp q 0 u t p f 0 u, rc l ¡ c r s tp q 0 u t0u .

(1.32)

• if δ $ 0 and δ $ N s {2, tp q δ u is complex and satisfies

Z 0 0 0 Z 0 & 4 Re pp q δ q Im pp q δ q B 5 Rep p f δ q Imp p f δ q C , c l ¡ c r cospδ αq c r sinpδ αq ¡c r sinpδ αq c l ¡ c r cospδ αq & 4
Re pp q δ q Im pp q δ q B 4 0 0 B .

(1.33)

• if δ N s {2, tp q Ns 2 u is real and satisfies

Z 0 $ 3 p q Ns 2 A t p f Ns 2 u, rc l c r s 3 p q Ns 2 A t0u . (1.34)
When all the individual responses have been computed, the full response is obtained by Fourier synthesis thanks to Eqn. (1.9).

The eigenmodes of the tuned disk are computed by imposing a zero forcing. It is equivalent to impose t p f δ u t0u for any value of δ. A consequence of this is that the eigenmodes of a tuned disk are strictly mono-harmonic and are fully determined by the resolution of the individual problems described in Eqns. (1.32), (1.33) and (1.34) without the right-hand terms. The modeshapes are then restored separately on the full disk according to the harmonic coefficient by

tqu rE 0 I N 0 s tp q 0 u if δ 0, tqu E Ns 2 I N 0 % 3 p q Ns 2 A if δ N s {2, tqu rE δ I N 0 s 4 Repp q δ q Impp q δ q B otherwise. (1.35)
It is classical to plot the eigenfrequencies of a tuned disk with respect to the Fourier harmonic coefficient δ, as it is done in Fig. 1.5 for a sample disk with 23 sectors whose geometry is similar to that of Fig. (i.) modes that correspond to disk-dominated motion, for which the sector to sector coupling ratio is high enough so that the strain energy spreads over the blades and the disk. Many papers cite curve veering [Yang and Griffin 2001, Judge et al. 2001, Bladh et al. 2002, Baik et al. 2004, Kenyon et al. 2005] to describe the behaviour of the frequencies associated with two different families of modes that tend to come close to each other. However, frequency veering is associated with frequencies that vary with respect to a continuous parameter [START_REF] Balmès | High modal density, curve veering, localization: a different perspective on the structural response[END_REF]], which is not the case for δ. Nevertheless, it is possible to use continuation techniques to zoom into some regions of assumed veering and to compute the radius of curvature of the continuous lines in the veering region to give an indication of blade-to-disk coupling ratio [Bladh et al. 2002, Baik et al. 2004].

(ii.) modes that correspond to blade-dominated motion, for which the sector to sector coupling ratio is relatively low, which tends to confine the strain energy to the blades [Ewins 1973, Bladh et al. 2002, Castanier and Pierre 2006]. This phenomenon can be seen as the transformation of backward-and forward-travelling waves reflecting at the inter-sector interfaces into standing waves [Touratier 1986,Castanier andPierre 2006]. Eigenfrequencies are distributed over horizontal lines whose asymptotes are given by the eigenvalues of the initial sector S 0 when its right and left interfaces I 0 r and I 0 l are assumed fixed. As a result, in a narrow frequency band whose upper limit is a "blade frequency" (more precisely, the frequency of such a fixed interface mode), there may be numerous modes corresponding to moderate to high Fourier harmonic coefficients. These regions where the modal density is high are referred to as accumulation bands, thus the frequencies of fixed interfaces modes are also often called accumulation frequencies. The families of modes are named after the blade modeshape associated with the accumulation with the pattern αTyβ where Ty describes the modeshape (F: flapwise bending, T: torsion, M: in-plane bending, E extension-compression and S: stripwise bending), α is the order of appearance in the family type which is the number of the main nodal lines (along the blade leading or trailing edge for F and along the blade tip edge for T and S, or orthogonal to the blade mean plane for M) and β is the number of the secondary nodal lines (along the blade leading or trailing edge, for S only). Note that all these shapes may not be observable depending on the geometry of the blade.

Such a figure provides very interesting information about why this decomposition presents a very significant gain in terms of numerical resolution. Indeed, it can be seen that the modal density is very high in certain frequency bands, in particular in the vicinity of the accumulation frequencies. As a result, considering one Fourier coefficient δ at a time leads to frequencies that are well separated (although they come in pairs for δ $ 0 and δ $ N s {2) which alleviates the convergence problems in the subspace iteration algorithms (e.g. Lanczos) used to solve the full finite element problem [Yang and Griffin 1997].

When it comes to compute the response to an external excitation, the procedure consists in a decomposition of the load into its Fourier harmonics, the resolution of as many subproblems as there exist such harmonics and the synthesis of the response from the harmonics of the displacement. For a mono-harmonic excitation, only the corresponding harmonic in the displacement will contribute to the response as displayed in Fig. 1.6 where a stationary external load with δ 2 is applied to a stationary bladed disk. Note that the two peaks are indistinguishable.

Remark: engine orders

The harmonic coefficients exhibited by the Fourier transform of the forcing excitation are commonly called "engine orders". These excitations are generally stationary, whereas the disk is rotating. In the rotating frame, the disk sees them as backward-travelling waves. If a mono-harmonic external load is applied with a harmonic δ, a given point on the disk will experience δ times the same forcing function per revolution of the disk, hence an excitation frequency of ω δ Ω{p2 πq (if the rotation speed Ω is expressed in rad.s ¡1 ).

Specific cases are:

• inertia (δ 0), • gravity (δ 1),

• single harmonic excitations due to the steady flow passing the upstream nozzle guide vanes (so-called engine order excitation with a unique δ),

• single point excitations (Dirac's comb).

Constraint elimination

When a constraint equation similar to Eqn. (1.22) is present in parallel to a dynamical equation like Eqn. (1.29), it states that some DOF are either imposed (single point constraints) or redundant (multi-point constraints). In this case, there is a need for DOF elimination. To do so, one builds the kinematically admissible subspace which is the kernel of the operator whose rcs is the matrix representation. This leads to a set of vectors rT Ker s that spans this kernel. The problem to be solved is then (1.37)

In the case of a tuned disk and due to the particular forms of matrices rZs and rcs (hence rT Ker s), the elimination process does not alter the particular form of the projected matrix rE t I Ñ 0 s r Zs rE I Ñ 0 s, where Ñ 0 is the number of unconstrained DOF in sector S 0 . If one defines r p T δ s as the vectors than span the kernel of rp c δ s present in Eqn. (1.19), the matrix with eliminated DOF is projected such that

E I Ñ 0 $ t r Zs E I Ñ 0 $ " " " " " " " " " ! p T t 0 Z 0 p T 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ p T t δ Z 0 0 0 Z 0 & p T δ ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ p T t Ns 2 Z 0 p T Ns 2 ( 0 0 0 0 0 0 0 0 0 )
.

(1.38)

A conclusion to be drawn from this equation is that if the constraint equation (1.22) is spatially periodic, the elimination process does not combine Fourier harmonics. Note that Eqn. (1.38) could have been obtained by performing separately the elimination of the constraints related to each harmonic in the projected matrices r p Z δ s in Eqns. (1.31).

This result has many important implications when a more global constraint equation is considered, especially when it includes not only the inter-sector continuity but also other conditions such as single point constraints or multi-point constraints that can be non-cyclic.

The lumped model described by Fig. 1.2 is used as an illustration. After the assembly of the sector matrices given by Eqn. (1.13) and after the elimination of the redundant left and right interface DOF, the resulting mass and stiffness matrices r M s and r Ks are projected onto

E I Ñ 0 $
and plotted in Fig. 1.7. These matrices are indeed block-diagonal and their shape is exactly that described by Eqn. (1.38), even if the diagonal blocks are no longer identical, as they were in Fig. 1.3. In structural dynamics, some norms are of common use (for details, see Appendix C). If rAs is a square matrix that is symmetric, positive and definite (resp. semi-definite), an inner product (resp. a semi-inner product) denoted x. , .y A can be defined such that dptqu, tpuq R N ¢ R N , xq , py A tqu t rAs tpu .

( 1.39) This allows to define a norm (resp. a semi-norm) denoted ||.|| A , which is given by dtqu R N , ||q|| A ptqu t rAs tquq 1 2 .

(1.40) Some choices for rAs exhibit the most commonly used norms and (semi-)inner products: (i.) if rAs rMs, since the mass matrix is definite, x. , .y M is an inner product and ||.|| M is a norm, called the norm in kinetic energy. The amount of kinetic energy that lies in the structure is given by dtqu R N , W M 1 2 ||q|| 2 M 1 2 tqu t rMs tqu .

(1.41) (ii.) if rAs rKs, since the stiffness matrix is semi-definite, x. , .y K is a semi-inner product and ||.|| K is a semi-norm, called the semi-norm in strain energy. The fact that there is no separation condition for a semi-norm, in other words that ||q|| K 0 does not necessarily mean that tqu t0u, has a physical meaning here: since any DOF vector describing a rigid body motion belongs the kernel of rKs, its norm in strain energy is zero, whereas it is obvious that it corresponds to a non-zero displacement. The amount of strain energy that lies in the structure is given by (iv.) a slightly different norm is also employed. The infinity-norm of R N is defined such that dtqu R N , ||q|| V maxptquq.

dtqu R N , W K 1 2 ||q||
(1.43)

Particular attention must be paid to the fact that all components must be of the same nature (translation or rotations) to allow comparisons in amplitude. Thus, it is mostly used in generalized modal analysis, where generalized DOF are often similar quantities, and in algorithms of optimal sensor placement with a single class of sensors [START_REF] Balmès | Orthogonal Maximum Sequence Sensor Placements Algorithms for modal tests, expansion and visibility[END_REF]].

Equation. (B.22) provides an interesting way to explicit the inner products and norm defined in Eqns. (1.39) and (1.40) in the case of tuned cyclic symmetric structures, that is to say where rAs is derived from rA 0 s such that rAs

I Ns A 0 $
. pxRe pp q δ q , Re pp p δ qy A 0 xIm pp q δ q , Im pp p δ qy A 0 q e p q Ns 2 , p p Ns

2 i A 0 .
(1.45)

The norm of any vector tqu is consequently

||q|| A ||p q 0 || A 0 Ns 2 ¡1 δ1 p||Re pp q δ q|| A 0 ||Im pp q δ q|| A 0 q § § § § § §p q Ns 2 § § § § § § A 0 .
(1.46)

Equation (1.46) is nothing else but Parseval's theorem relative to the inner product x. , .y A (see Appendix B. 1.3). Recall that rEs is a unitary matrix, therefore it preserves the inner products, hence the norms.

This result is fundamental since the total energy can be recovered by summing the energy associated with each Fourier harmonic, with a coefficient equal to 1 thanks to the unitary Fourier transform but with distinct contributions of the real and imaginary parts. Conversely, the fraction of strain energy associated with each Fourier harmonic can be computed.

Another consequence of this is a new definition of the Modal Assurance Criterion (MAC) to check the orthogonality of complex mono-harmonic vectors after they have been computed by the means of Eqns. (1.32) to (1.34). If tp q δ u and tp p δ u are two such vectors, this definition is MAC A 0 pp q δ , p p δ q |xRe pp q δ q , Re pp p δ qy A 0 xIm pp q δ q , Im pp p δ qy A 0 | 2 p||Re pp q δ q|| A 0 ||Im pp q δ q|| A 0 q p||Re pp p δ q|| A 0 ||Im pp p δ q|| A 0 q .

(1.47)

This updated definition enables an auto-MAC where each vector in a pair of complex conjugate shapes is matched with itself only and not also to the other vector in the pair.

Comments on the methodological aspects

This section can be summarized by the flowchart in Fig. 1.8.

Mistuning

Sec. 1.1 describes the case of a geometrically periodic structure whose mechanical properties are identical from one bladed sector to another. Such a property is very convenient since it allows to decompose the finite element problem expressed for the whole structure into independent sub-problems driven by the Fourier harmonic coefficient δ. The question is then to know what happens when this assumption fails, in other words when the disk is said to be mistuned.

Characteristics of mistuning

Indeed, any geometrically periodic structure is said to be mistuned when there exists discrepancies in the geometry or in the mechanical properties from a bladed sector to another. In fact, two different words exist, detuning and mistuning, related to the fact that this phenomenon is intentional or unintentional, respectively. Therefore, detuning can result in manageable patterns of varying parameters (sines, alternations of two types of blades...) whereas mistuning is generally considered as a random phenomenon due to its origins (that remain not well known). Numerous studies on this subject exist in the literature, only the most significant ones are cited here.

A set of parameters whose spatial variations have a significant impact on the behaviour of the nominally tuned structure was identified from various lumped models similar to those shown in Fig. 1.9 as well as simple or realistic disk models. Figure 1.9: Lumped model of a single mistuned disk These identified parameters are:

• the mass of the structure [Judge et al. 2001, Mignolet et al. 2001, Rivas-Guerra et al. 2001, Yang and Griffin 2001, Lalanne 2005, Petrov and Ewins 2005, Christensen and Santos 2005],

• the stiffness of the disk or of the blades [Mignolet and Lin 1997, Yang and Griffin 1997, Whitehead 1998, Mignolet et al. 2000a, Mignolet et al. 2000b, Judge et al. 2001, Yang and Griffin 2001, Mignolet et al. 2001, Rivas-Guerra et al. 2001, Castanier and Pierre 2002, Moyroud et al. 2002, Bladh et al. 2003, Baik et al. 2004, Ratier 2004, Lalanne 2005, Petrov and Ewins 2005],

• the damping [Lin and Mignolet 1996, Mignolet and Lin 1997, Petrov and Ewins 2005],

• the whole dynamic stiffness [Petrov et al. 2002, Rivas-Guerra and[START_REF] Rivas-Guerra | [END_REF]],

• the natural frequencies of the cantilevered blades [Whitehead 1966, Ewins 1969, Griffin and Hoosac 1984, Bladh et al. 2001a, Bladh et al. 2001b, Bladh et al. 2002, Pierre and Murthy 2002] or of a bladed sector whose inter-sector interfaces are fixed [Feiner and Griffin 2002, Feiner and Griffin 2004a, Feiner and Griffin 2004b], which can be interpreted as a combined effect of the mass and the stiffness (of the blades in the first case and of the bladed sector in the second),

• the aeroelastic coupling [Lalanne andTouratier 2000, Pierre and[START_REF] Pierre | [END_REF],

• the geometry [Ewins 1973, Ewins 1976], which can originate from the manufacturing process [Feiner andGriffin 2004b, Capiez-Lernout et al. 2005], from in-service wear or from an accidental scenario (a foreign object that impacts the blade which at best leads to a permanent deformation of the blade and at worst to a blade-off event and further damage to the downstream stages) [Ewins 1973, Lim et al. 2004].

When the variations in the mechanical parameters or in the geometry are small enough to neglect the changes in the cantilevered blade modeshapes, mistuning is said to be small and large otherwise. In the latter case, the assumption of geometrical periodicity fails and it will not be considered in the present chapter, even if some papers allow to deal with a large mistuned structure with the help from the underlying cyclically symmetric system [Rixen and Lohman 2005].

Whatever its amplitude, any loss of symmetry induces a split of the pairs of complex conjugate modes. A zoom on the peaks of the forced response function of a mistuned disk to a single engine order excitation reveals two peaks (see Fig. 1.10a), all the more separated that the damping ratio is small [Ewins 1973] which was demonstrated experimentally [Ewins 1976, Judge et al. 2001]. All the mistuned frequencies arise from a scattering of the tuned frequencies [Ewins 1969, Ewins 1973, Ewins 1976, Griffin and Hoosac 1984]. There arises the numerical issue associated with the convergence of the subspace iteration algorithm towards near-pairs of frequencies that is bypassed in the tuned problem by the help from the decomposition into independent sub-problems.

Many published studies reveal that mistuning may affect the amplitude in the response of the blades and subsequently the stress levels they encounter. The amplification factor varies according to the underlying model [Whitehead 1998,Mignolet et al. 2000a,Mignolet et al. 2000b,Petrov and Ewins 2005, Han and Mignolet 2007, Han et al. 2007]. Besides, the modal controllabilities and observabilities strongly depend on the mistuning ratio [Christensen and Santos 2005]. Therefore, it is mandatory to design an indicator that would be more accurate than simply comparing amplitudes of the displacement according to mistuning, for instance by considering stress levels in the blades [START_REF] Lim | [END_REF]]. However, mistuning is not only a phenomenon with adverse effects since it provides structures with stabilization towards flutter [Srinivasan 1997, Lalanne and[START_REF] Lalanne | [END_REF].

Mistuning, in opposition to detuning, is by essence a random phenomenon [START_REF] Griffin | [END_REF]. A first approach consists in choosing a few parameters, generating Monte-Carlo samples and computing the response level for each sample. The number of samples can be very high to accurately account for the extremities of the distribution laws. Reduced statistical laws such as that of Weibull's can be used [Zaretsky et al. 2003] to decrease that number of samples. Such methods require a knowledge of the probability densities of the involved parameters. Conversely, from experimental data collected for a set of nearly identical structures, it is possible to return to the characteristics of the nominally tuned structure [Mignolet and Lin 1997] as well as those of each mistuned structure [Feiner andGriffin 2004a,Feiner andGriffin 2004b]. To eliminate the dependence on explicit parameters, a non-parametric approach was proposed to take into account a more global uncertainty, of the parameters as well as of the model itself, by considering random variations of the mean finite element matrices and making a statistical analysis of the response levels [START_REF] Capiez-Lernout | Dynamique des structures tournantes à symétrie cyclique en présence d'incertitudes aléatoires[END_REF][START_REF] Capiez-Lernout | Nonparametric Modeling of Random Uncertainties for Dynamic Response of Mistuned Bladed Disks[END_REF], Capiez-Lernout et al. 2005].

In the tuned case, the energy spreads over all blades. When the disk is mistuned, the response (hence the energy) may be localized to a subset of blades, even a single one [Whitehead 1966,Pierre and[START_REF] Pierre | [END_REF]. Peak scattering occurs, as depicted in Fig. 1.10 where the disk of Sec. 1.1 has its 23 blades mistuned in frequency with a standard deviation of 0.5%, 5%, 10% and 20% and is subject to the same δ 2 excitation. It is essentially degenerated tuned modes that are observed in that very case. Localization strongly depends on the coupling between the different parts of the structure (blade-to-blade, disk-to-blade) [Whitehead 1966, Ewins 1973, Xie and Ariaratnam 1996a, Xie and Ariaratnam 1996b, Mignolet et al. 2000a, Mignolet et al. 2000b, Baik et al. 2004, Castanier and Pierre 2006]. Indeed, the localization ratio increases as the blade-to-blade coupling ratio decreases: the strain energy has difficulty spreading over the adjacent blades (which must not be considered as local damping) [Castanier and Pierre 2006]. To study this phenomenon, partially mistuned systems were considered so that a single mistuned blade vibrates and is observed from its neighbours [Rivas-Guerra and [START_REF] Rivas-Guerra | [END_REF]. Conversely, if the coupling is strong the strain energy spreads easily over the full disk and the mistuned behaviour is close to the tuned one [Castanier and Pierre 2006]. The disk-to-blade coupling also drives the modal density in the regions of disk-dominated motion, therefore it is necessary to look carefully at these regions [Ewins 1973[START_REF] Balmès | High modal density, curve veering, localization: a different perspective on the structural response[END_REF], Bladh et al. 2002, Lalanne 2005].

Since mistuning is hardly manageable, it was proposed to artificially detune a bladed disk with either partial or alternating patterns and look at how the detuned system is sensitive to mistuning [START_REF] Castanier | Using Intentional Mistuning in the Design of Turbomachinery Rotors[END_REF]]. An optimization technique based on a genetic algorithm was proposed to look for the best pattern given only two different types of blades [Choi et al. 2003, Han and[START_REF] Han | [END_REF]. Note that the pattern is only optimal for a single loading condition, which limits its practical use.

Note finally that all these studies consider the disk at rest. The inertial loading tends to stiffen the disk (see Chap. 3) which tends to limit the effects of mistuning [Moyroud et al. 2002].

Phenomenological aspects

As described in the previous section, mistuning is either a physical phenomenon, related to variations in the mechanical parameters from a bladed sector to another, or a geometrical phenomenon when the nominal geometry of one or more blades are altered. The second case is outside the framework of this chapter since changes in the geometry implies that there is no longer a series of corresponding points px s q s with an occurrence on each blade, that is to say the hypothesis of cyclic symmetry no longer holds (see Sec. 1.1.1). It can however be treated by subsequent reduction techniques based on substructuring procedures (see Chap. 2). However, in the first case, since the geometry is left unchanged, mistuning only affects the physical parameters, hence the finite element matrices.

Let rAs be a tuned matrix that depends on one parameter. From Sec. 1.1.3, this matrix depends only on the matrix rA 0 s of sector S 0 and is rAs

I Ns A 0 $
, before any DOF elimination. Variations in parameters from a bladed sector to another are assumed fully contained in a matrix r∆As that usually has no characteristic shape except that it contains one block per sector denoted r∆A s s, with at least one non zero block. According to the mistuning pattern applied to the parameter, these blocks can be associated with only a few blades (partial mistuning), or if all sectors are concerned, they can be randomly distributed (random mistuning) or harmonically distributed (harmonic mistuning) over the diagonal. A more explicit form of rA ∆As is given by rA ∆As .48) The projection of matrix rA ∆As onto rE

" " " " " " ! A 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ A 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ A 0 ( 0 0 0 0 0 0 ) " " " " " " ! ∆A 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ ∆A s ¤ ¤ ¤ 0 . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ ∆A Ns¡1 ( 0 0 0 0 0 0 ) . ( 1
I N 0 s is rE I N 0 s t rA ∆As rE I N 0 s r p A ∆ p As, (1.49) r p
As is given by Eqn. (1.12) and still do not couple harmonics together. The problem is that r∆ p As is not a block-diagonal matrix and may couple many, if not all, harmonics.

An interesting case is when the matrix r∆As derives from a harmonic mistuning pattern. In other words, the diagonal blocks r∆A s s are distributed according to a sine/cosine pattern driven by ζ similarly to the Fourier decomposition of the displacement (with ζ used here so that there is no possible confusion with the Fourier harmonic coefficients δ that may all exist) and each block r∆A s s is ds v0 Alternate mistuning is then achieved with the cosine pattern when ζ N s {2, in other words two adjacent sectors are mistuned with the same matrix r∆AN s 2 s but with opposite signs. Random mistuning is a more general case which involves all possible sine/cosine patterns with ζ in v0, N s {2w summed to form rAs. As a result, the conditions described in Eqn. (1.53) are fulfilled whatever the pair pδ, δ I q so that all harmonics are coupled together.

δ δ I ζ rN s s, N s ¡ δ δ I ζ rN s s, δ N s ¡ δ I ζ rN s s, 2 N s ¡ δ ¡ δ I ζ rN s s.
An illustration of this is obtained if one modifies the lumped model in Fig. 1.2 by letting the masses and stiffnesses vary from sector to sector, as described in Fig. 1.9.

Figure 1.11 is obtained for a disk with 12 sectors and 4 springs per sector. In each case, the matrices of the full disk are assembled first and DOF elimination is then performed. Figure 1.11a is presented as a reminder of the tuned case. In Fig. 1.11b, a sine harmonic mistuning pattern with ζ 1 is applied to the stiffness matrix so that k s d k 0 d p1 .05 sinps ζ αqq. This figure shows that only the pairs of harmonics pδ, δ I q such that δ δ I 1 (modulo the aliasing) are coupled together, which corresponds to the blocks located on first upper and lower diagonals, and that others harmonics remain uncoupled. The blocks on the diagonal correspond to r p As. In Fig. 1.11c, a cosine harmonic mistuning pattern with ζ 6 is applied to the stiffness matrix so that k s d k 0 d p1 p¡1q s .05q, which corresponds to alternate mistuning. In this case, only the pairs of harmonics pδ, δ I q such that δ δ I 6 (modulo the aliasing) are coupled together, which corresponds to the blocks located on the anti-diagonal. It implies that there is a block on the diagonal with δ δ I N s {4. Once more, other harmonics remain uncoupled. Finally, in Fig. 1.11d, a random mistuning pattern where parameters in r∆Ks vary between ¡5% and 5% of the mean value is added to the tuned stiffness matrix. It can clearly be seen that all harmonics are coupled together since all sine and cosine patterns are involved.

These results point out the need for full finite element computations in the case of matrix mistuning, since mistuning gives rise to intra-disk harmonic coupling. The mistuning pattern drives the number of harmonics involved. Among those, alternate mistuning is a particular case where Fourier harmonics of the discrete displacement are coupled by pairs without overlapping as can be seen in Fig. 1.11c. Therefore, separate problems involving two harmonics could be solved. Such a case is rarely seen in practice where most of the cases deal with random mistuning. Full disk strategies, but not necessarily involving the full 3D mesh of the structure, have to be designed.

Multi-stage dynamics

Sections 1.1 and 1.2 consider a single rotationally periodic structures with various patterns for the distribution of its mechanical properties. The question is now to determine the behaviour of assemblies of such cyclically symmetric structures so that full rotor models can be examined.

In industry, different methods already exist according to the expected degree of accuracy of the models. They are presented here in their order of appearance: (i.) the simplest approach consists in assuming that the rotor is axisymmetric. A longitudinal section is meshed with axisymmetric Fourier 2D elements as displayed in Fig. 1.12a. The main drawback of this method is that the geometry of the blades is totally bypassed. Such a method is only suited to rotors with a strongly disk-dominated motion, which is rarely the case in aero engines.

(ii.) the second one focuses on a given disk inside its cyclic symmetry group (see Sec. 1.1).

Adjacent stages are added so that they virtually belong to this symmetry group by modifying their blade distribution and even simplifying their geometry, as depicted in Fig. 1.12b. This method provides good results for the concerned stage in terms of coupling but has to be repeated until all stages have been analyzed.

p n ¹ p P p ¤ n! ¢ 1 ¡ 1 p 2
, where P N is the set of the prime numbers. Cesàro demonstrated that such a sequence has a limit of 6{π 2 when n increases towards V. The probability that the number of blades of two adjacent disks are coprime is 60.88 % under 5! 120 blades and stays beyond 6{π 2 60.79 % whatever these two numbers. Since the number of blades in each stage is a constraint induced by aerodynamic efficiency, such convenient cases are exceptional in practice. A pair of even numbers is the most frequent case, it allows to consider one half of the rotor.

(iv.) when the latter solution is not possible, a simplified model of the whole geometry using both 3D solid and 2D shell elements can be used, provided that such simplifications do not affect the behaviour too much, in particular at the connections between shells and solid elements. Such a model is presented in Fig. 1.12d.

(v.) more recently, the classical approach of cyclic symmetry for single bladed disks was extended to multi-stage assemblies so that the eigensolutions associated with a single Fourier harmonic coefficient δ can be computed under specific assumptions [Laxalde et al. 2006, Laxalde et al. 2007a]. A typical solution is presented in Fig. 1.12e.

The Multi-Stage Cyclic Symmetry approach is the core of this section. The variant described here derives from the original formulation which considers joint disks with or without conforming meshes. It is adapted to the particular framework in which the cyclically symmetric disks are connected with a volumic ring which concentrates the asymmetry.

Assembly of two disks through an intermediate ring

Let us consider two bladed disks D 1 and D 2 that are connected through a ring R of solid elements such as presented in Fig. 1.13. Disk D 1 and ring R share the same surface I 1 a . Ring R and disk D 2 share the same surface I 2 f . Due to the way it is meshed, the inter-disk ring R has a single element along the z direction, thus it has no interior DOF when it is made of first order elements. D 1 and D 2 have N 1 and N 2 DOF respectively. These two DOF sets are disjoint. The

total number of DOF is N N 1 N 2 . E c E e x e y e z f f f f f x ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ f f f w f f f f f x ¢ ¢ ¢ ¢ ¢ ¢ D 1 R D 2 I 1 f I 1 a I R f I R a I 2 f I 2 a
R u t q R f q R a u t .
Disks D 1 and D 2 are geometrically periodic. Therefore, from Sec. 1.1, Fourier harmonics can be defined separately for each disk thanks to Eqn. (1.9) such that for d t1, 2u, tp q d u r p E d I N d,0 s tq d u.

(1.54)

N d,0 is the number of DOF (without elimination) that belong to sector S d,0 . The discrete displacement is recovered by means of Eqn. (1.9) such that for d t1, 2u, tq d u rE d I N d,0 s tp q d u.

(1.55)

Any matrix rAs relative to this rotor comes from the assembly of matrices rA 1 s, rA 2 s and rA R s defined on D 1 , D 2 and R respectively such that rAs

A 1 0 0 A 2 & A R $ . (1.56)
Matrix rAs is projected onto the subspace spanned by the columns of the block-diagonal matrix

rdiagpE 1 I N 1,0 , E 2 I N 2,0 qs r p As E 1 I N 1,0 0 0 E 2 I N 2,0 & t rAs E 1 I N 1,0 0 0 E 2 I N 2,0 & . (1.57) More precisely, r p As is r p As p A 1 0 0 p A 2 ' r p A R s, (1.58) with r p A R s E 1 I N 1,0 0 0 E 2 I N 2,0 & t A R $ E 1 I N 1,0 0 0 E 2 I N 2,0 & .
(1.59)

The topology of r p A d s, d t1, 2u, depends on the degree of mistuning of the disk D d , as was shown in Secs. 1.1.3 and 1.2.2. In the particular case when disk D d is tuned, matrix rA d s is given by rI Ns d A d s. Consequently, matrix r p A d s is given by Eqn. (1.12) and it does not couple the Fourier harmonics tp q d δ u together.

Matrix r p

A R s is somewhat more complex since ring R does not satisfy cyclic symmetry. If this matrix is decomposed into blocks indexed by dd I that couple disks D d and D d 1 , these blocks are

r p A R dd 1 s rE d I N d,0 s t rA R dd 1 s rE d 1 I N d 1 ,0 s. (1.60)
These blocks can themselves be subdivided into blocks indexed by δδ I that couple Fourier

harmonic tq d δ u of D d to Fourier harmonic tq d 1 δ 1 u of D d 1 , leading to r p A R dd 1 ,δδ 1 s rE d δ I N d,0 s t rA R dd 1 s rE d 1 δ 1 I N d 1 ,0 s. (1.61)
Equation (1.61) shows that:

(i.) the Fourier harmonics of a single disk are not coupled together, therefore the ring matrix does not contribute to intra-disk harmonic coupling and only mistuning may give rise to it, (ii.) the Fourier harmonics of a given disk are coupled to those of the other disk and this is the only contribution to inter-disk harmonic coupling. The lumped model used for illustrations in the mono-disk sections is enriched to cover multi-stage cases. To this extent, two tuned lumped models are coupled together by means of springs, as displayed in Fig. 1.14. For the sake of simplicity, it is assumed that disk D 1 (resp. disk D 2 ) has N 2 s (resp. N 1 s ) springs per sector. After assembly of mass and stiffness matrices of the rotor and DOF elimination, matrices r x M s and r p Ks are computed. Figure 1.15 concerns two lumped disks with 4 and 8 sectors respectively. It can be seen in Figs. 1.15a and 1.15b that the matrices of the individual disks are indeed block-diagonal. Notice that harmonics are grouped together with the coloring convention of Fig. 1.3 and disks are interlaced so that the row/column pattern is pδ 0, D 1 q, pδ 0, D 2 q, pδ 1, D 1 q, pδ 1, D 2 q... Blocks of disk D d are not present as soon as δ ¡ N d s {2. Figure 1.15c however shows aliasing in the off-diagonal coupling terms which leads to inter-disk harmonic coupling. δ 3 is indeed δ 3 on disk D 2 but δ 4 ¡ 3 1 on disk D 1 , thus one sees non-zero coupling in the ptp q 1 1 u, tp q 2 3 uq blocks. δ 4 is indeed δ 4 on disk D 2 but δ 4 ¡ 4 0 on disk D 1 and one sees non-zero coupling in the ptp q 1 0 u, tp q 2 4 uq blocks. shows that there is no contribution to intra-disk harmonic coupling either. But inter-harmonic coupling involves many more harmonics than in the previous case. Indeed, all harmonics are coupled from a disk to another.

Hypothesis of decoupled harmonics

The inter-harmonic coupling shown in Eqn. (1.61) prevents any rotor problem from being decomposed into individual sub-problems for every harmonic coefficient as it could be done in mono-disk analyses, see Sec. 1.1.4. This also implies that free vibration modes of a rotor can only be mono-harmonic if the inter-harmonic coupling coefficients in (1.61) have a negligible influence.

Since decoupled harmonics present a very high interest in practice, an approximation can be made. In the following, one assumes that in any matrix r p As as defined in Eqn. (1.56):

• the sub-matrices indexed by δδ I in the blocks r p A R dd s are neglected if δ $ δ I , • the sub-matrices indexed by δδ I in the blocks r p A R dd 1 s, d $ d I , are neglected except the pairs pδ, δ I q that satisfy any of the following conditions

δ δ I rN s m s, δ N s m ¡ δ I rN s m s, (1.62)
where

N s m minpN d s , N d 1 s q. (1.63) With this approximation, r p A R δ s is defined retroactively such that r p A R δ s E 1 δ I N 1,0 0 0 E 2 δ I N 2,0 & t rA R s E 1 δ I N 1,0 0 0 E 2 δ I N 2,0 & . (1.64)
Theoretically, the definition of Eqn. (1.64) is valid for any Fourier harmonic coefficient δ even if it is higher than N d s {2 for any disk. In this case, subscript δ appears in each matrix rE d δ s

modulo the aliasing with a periodicity of N d s with the correspondence

if k N d s δ p2k 1q N d s {2 , k N , rE d δ s rE d δ¡k N d s s, if p2k 1q N d s {2 δ pk 1q N d s , k N , rE d δ s rE d pk 1q N d s ¡δ s.
(1.65) In the case of two disks, this correspondence phenomenon occurs with a global periodicity of lcmpN d s , N d 1 s q (least common multiple of N d s and N d 1 s

). An example where disk D 1 has 4 or 5 sectors, and disk D 2 has 8 sectors is displayed in Fig. 1.17. In this figure, the two curves start from δ 0 and meet at δ 8 lcmp4, 8q and δ 40 lcmp5, 8q respectively. Besides, notice that these curves are symmetric with respect to δ 4 and δ 20 respectively, which implies that δ and lcmpN d s , N d 1 s q ¡ δ are similar. For more than two disks, the situation becomes more complex since more curves have to cross at a single point. This is illustrated by Fig. 1.18 relative to the lumped rotor models studied in Figs. 1.15 and 1.16.

In this figure, blocks of matrix r p K R s, given by Eqn. (1.59), are virtually identified to the blocks in r p K R δ s when δ varies from 0 to lcmpN 1 s , N 2 s q{2. In Fig. 1.18a, since lcmp4, 8q 8, only a few blocks are filled since the harmonic correspondence pattern is very simple. On the contrary, in

Response of a rotor made of tuned disks to external excitations

A general dynamic problem is formulated for the rotor rZs tqu tfu , (1.66) with tqu 2 q 1 q 2 @ and tfu

2 f 1 f 2 @
. It can be rewritten r p Zs tp qu t p f u.

(1.67)

This equation is valid whatever the degree of mistuning of each disk. Recall that the latter only drives the degree of intra-disk coupling between the Fourier harmonics. In Sec. 1.3.1, it was demonstrated that in the general case, harmonics are coupled together even though all disks are tuned, which precludes a decomposition of the mechanical problem into mono-harmonic sub-problems. Nevertheless, such a decomposition is achieved under the assumption of decoupled harmonics described in Sec. 1.3.2. In the case of a rotor made of tuned disks, Eqn. (1.67) results in a family of parametrized problems whose solutions are spatially mono-harmonic:

• δ 0 rlcmpN 1 s , N 2
s qs. The restrictions of the mono-harmonic solution to disks D 1 and D 2 are both real. Such solutions satisfy

¢ Z 1,0 0 0 Z 2,0 & p Z R 0 % 4 p q 1 0 p q 2 0 B 5 p f 1 0 p f 2 0 C , c 1 l ¡ c 1 r 0 0 c 2 l ¡ c 2 r & 4 p q 1 0 p q 2 0 B t0u .
(1.68)

• δ % 0 rlcmpN 1 s , N 2
s qs and δ % 0 rN 1 s {2s and δ % 0 rN 2 s {2s. The restrictions of the monoharmonic solution to disks D 1 and D 2 are both complex. Such solutions satisfy

¤ ¦ ¦ ¥ " " ! Z 1,0 0 0 0 0 Z 1,0 0 0 0 0 Z 2,0 0 0 0 0 Z 2,0 ( 0 0 ) p Z R δ % 6 9 9 8 9 9 7 Repp q 1 δ q Impp q 1 δ q Repp q 2 δ q Impp q 2 δ q D G G F G G E 6 9 9 9 8 9 9 9 7 Rep p f 1 δ q Imp p f 1 δ q Rep p f 2 δ q Imp p f 2 δ q D G G G F G G G E , " " ! c 1 l ¡ cospδ α 1 q c 1 r sinpδ α 1 q c 1 r 0 0 ¡ sinpδ α 1 q c 1 r c 1 l ¡ cospδ α 1 q c 1 r 0 0 0 0 c 2 l ¡ cospδ α 2 q c 2 r sinpδ α 2 q c 2 r 0 0 ¡ sinpδ α 2 q c 2 r c 2 l ¡ cospδ α 2 q c 2 r ( 0 0 ) 6 9 9 8 9 9 7 Repp q 1 δ q Impp q 1 δ q Repp q 2 δ q Impp q 2 δ q D G G F G G E t0u .
(1.69)

• δ % 0 rlcmpN 1 s , N 2
s qs and δ 0 rN 1 s {2s. The restriction of the mono-harmonic solution to

disk D 1 is real but its restriction to disk D 2 is complex. Such solutions satisfy ¤ ¥ ! Z 1,0 0 0 0 Z 2,0 0 0 0 Z 2,0 ( ) p Z R δ % 6 8 7 p q 1 δ Repp q 2 δ q Impp q 2 δ q D F E 6 9 8 9 7 p f 1 δ Rep p f 2 δ q Imp p f 2 δ q D G F G E , ! c 1 l c 1 r 0 0 0 c 2 l ¡ cospδ α 2 q c 2 r sinpδ α 2 q c 2 r 0 ¡ sinpδ α 2 q c 2 r c 2 l ¡ cospδ α 2 q c 2 r ( ) 6 8 7 p q 1 δ Repp q 2 δ q Impp q 2 δ q D F E t0u .
(1.70)

• δ % 0 rlcmpN 1 s , N 2
s qs and δ 0 rN 2 s {2s. The restriction of the mono-harmonic solution to

disk D 1 is complex but its restriction to disk D 2 is real. Such solutions satisfy ¤ ¥ ! Z 1,0 0 0 0 Z 1,0 0 0 0 Z 2,0 ( ) p Z R δ % 6 8 7 Repp q 1 δ q Impp q 1 δ q p q 2 δ D F E 6 9 8 9 7 Rep p f 1 δ q Imp p f 1 δ q p f 2 δ D G F G E , ! c 1 l ¡ cospδ α 1 q c 1 r sinpδ α 1 q c 1 r 0 ¡ sinpδ α 1 q c 1 r c 1 l ¡ cospδ α 1 q c 1 r 0 0 c 2 l c 2 r ( ) 6 8 7 Repp q 1 δ q Impp q 1 δ q p q 2 δ D F E t0u . (1.71)
The equations that come first in the problems above are not sufficient to define a mono-harmonic solution. The inter-sector continuity is added as a constraint equation that states how the shapes are dephased between two adjacent sectors, exactly like for the mono-disk case described in Sec. 1.1.5.

Mono-harmonic eigensolutions are derived from Eqns. (1.68) to (1.71) by letting the right hand terms be zero. One insists on the fact that these mono-harmonic eigensolutions may differ from the true normal modes of the rotor.

In the second problem, the projection of the dynamic stiffness matrix of any of the disks is equivalent to duplicate this matrix. This is unfortunately not the case for the ring, which leads to near-pair eigensolutions. For thinner meshes of the ring, the near-pair solutions would converge Like in the mono-disk case, the mono-harmonic eigensolutions can be gathered into a graph that displays them according to the value of the Fourier harmonic coefficient δ. This is done in Fig. 1.19 for an assembly of two disks with 23 and 45 blades respectively.

A few comments arise from the study of Fig. 1.19:

• Multi-stage effects are distinguishable provided that the eigenmodes of the individual disks are available. Figure 1.19b is built by considering each disk with either its inter-stage interface being free or fixed. As a result, this figure shows that inside a triangular area, some of the multi-stage eigensolutions are comprised between the fixed-and free-interface modes of one of the two disks which leads to think that, as expected, the real boundary condition lies somewhere between that two types. However, some other eigensolutions are not clearly related to any of the individual disk modes which emphasizes the need for true multi-stage computations.

• Accumulations are clearly seen in Fig. 1.19a. Since the modal density of a single disk may already be very high in the vicinity of its accumulation frequencies (see Sec. 1.1), this phenomenon is emphasized by the presence of all the accumulation bands of all stages.

More important is that most of the accumulations remain unaffected by multi-stage effects (especially at low frequencies), merely due to a small disk-to-blade coupling, which tends to confine the strain energy on the blades [Castanier and Pierre 2006]. consequently, the global motion is very close to that of the individual disks whatever their inter-stage boundary conditions since the latter do not affect the blade-dominated motion as can be seen in Fig. 1.19b. Accumulation with significant discrepancies are also observed, which leads to think that some modes accumulate towards modes that still have a significant amount of energy in the disk. This is coherent with the single-stage observations and are not multi-stage induced phenomena.

• Aliasing occurs for the disk with the smaller number of blades. In the present case where this number is 23, the mono-harmonic eigensolutions with δ ¡ 11 mostly localized to disk D 1 are aliases of solutions with δ 11. The question is to know whether such aliases are physical solutions or not. In particular, is it possible that an alias exists on both disks simultaneously? Can a high diameter solution exist both as part of the accumulation band of a given disk, thus potentially localized, and alias of a low-diameter solution, potentially coupled, of another disk? This issue is discussed in parallel to the study of the 360°reduced models.

Assembly of two disks with a constraint equation

In the previous sections, the two disks were supposed connected with an intermediate ring. (ii.) inter-disk continuity. It requires the definition of rc R s whose sparsity depends on the meshes of I 1 a and I 2 f . When meshes are not compatible, a given DOF of I 1 a (resp. I 2 f ) is the linear combination of DOF of I 2 f (resp. I 1 a ) and matrix rc R s involves projection operators.

The global constraint equation is then

rcs 4 q 1 q 2 B 4 0 0 B , (1.72) with rcs ! c 1 0 0 c 2 c R ( ) . (1.73)
It is demonstrated in Sec. 1.1.5 that, in the case of a tuned disk, the transformation of rc d s into rp c d s does not perturbe the tuning of the system, thus it does not contribute to intra-disk harmonic coupling either. However, the latter arises from the transformation of rc R s into rp c R s. This inter-harmonic coupling is propagated to the generating set of its kernel, namely r p T Ker s span pkerprp csqq, called the constrained subspace. In the case of the rotor with 4 and 8 sectors, the whole lumped model can be considered as a cyclically symmetric structure with 4 sectors. Consequently, three groups of independent columns are clearly seen and they generate the global harmonics, that is to say relative to the whole assembly, with δ t0, 1, 2u. The group of vectors in the middle (columns with δ 1) is twice longer than the others (columns with δ t0, 2u). Aliasing is therefore observed in the sense that:

• the global δ 0 is generated by δ 0 (disk D 1 ) and δ t0, 4u (disk D 2 ),

• the global δ 1 is generated by δ 1 (disk D 1 ) and δ t1, 3u (disk D 2 ), The lumped rotor with 5 and 8 blades does not belong to any particular symmetry group when considered as a whole. It is clearly seen that the constrained subspace is generated by vectors involving multiple harmonics. An assumption of decoupled harmonics can be made, however it has to be formulated relatively to Eqn. (1.72), in other words relatively to rcs.

• the global δ 2 is generated by δ 2 (disk D 1 ) and δ 2 (disk D 2 ) (a) (b)
Chapter 1 provides the means to deal with assemblies of bladed disks in various configurations. It makes it possible to compute (i.) all the eigenmodes of a single tuned disk by considering all possible Fourier harmonics, or a subset of them according to the spatial spectrum of the external excitation, (ii.) all the eigenmodes of a single mistuned disk by considering the full finite element problem because of inter-harmonic coupling, (iii.) all the mono-harmonic eigensolutions of a bladed rotor made of tuned disks by considering all possible harmonics even though (a) such solutions are not necessarily modes and (b) aliased solutions may exist, (iv.) all the eigenmodes of a bladed rotor made of tuned or mistuned disks by considering the full finite element problem .

The correspondence between modes and mono-harmonic eigensolutions holds only for an isolated tuned disk. Harmonic synthesis is possible in that case, that is to say, the response to an external excitation is the superposition of the individual responses to each single spatial harmonic in the load. Unfortunately, this capability no longer exists for mistuned disks nor multi-stage bladed rotors. Mistuning perturbs the overall symmetry of the system and modes are necessarily spatially multi-harmonic. When it comes to multi-stage assemblies of tuned disks, despite the convenience of the mono-harmonic eigensolutions, the presence of non-physical aliases among them prevents the synthesis of the response from a subset of such solutions, unless an accurate criterion to eliminate aliases is available, which is not as simple as it seems.

Mistuning first requires that distinct models of each individual sector have to be assembled into a 360°disk model (resp. rotor model). The main issue concerns the size of such models.

Typical meshes of a bladed sector have up to 10 5 ¡ 5.10 5 DOF due to the very high spatial resolution required to accurately reproduce the stress distributions in critical areas, e.g. the blades. Moreover, this issue is emphasized because such structures have a very high modal density, especially in the accumulation bands. The storage of the whole set of solutions would be close to the upper limit of the current hardware available. As an illustration, 1 000 modes for a structure with 15.10 6 DOF needs 200 Gb to be stored. In addition, when it comes to the statistical analysis of a family of mistuned disks or rotors, the single deterministic computation is replaced by a few thousand of Monte-Carlo samples. The need for compact but accurate reduced models is clear. This chapter aims to present two substructuring techniques designed to build reduced-order models of single-as well as multi-stage bladed rotors. These techniques are quite similar by nature in the sense that they associate a superelement to each bladed sector that aims to represent its dynamics with a controlled accuracy. These superelements are later assembled to form a 360°disk or rotor model.

Moreover, computing exactly the whole set of mono-harmonic eigensolutions of a tuned bladed disk (resp. a bladed rotor made of tuned disks) may represent a very high numerical cost in industrial applications where a stage can have up to 100 ¡ 150 blades. Storing such solutions, when disks have dense meshes as described above, would also require a prohibitive amount of disk space. Consequently, the two proposed reduction techniques rely on the concept of target vectors. The latter are mono-harmonic eigensolutions selected for a particular reason, for instance because the external load exhibits these harmonics in its spectrum. Vectors in this set are not linearly independent and worse, they can be very similar, especially in the accumulation bands. Orthogonalization and elimination procedures are specifically introduced to address this issue. The subsequent reduced model would therefore reproduce these solutions exactly (a) in the case of a single tuned disk and (b) in the case of a rotor made of tuned disks when it happens that a subset of modes is matched to a subset of mono-harmonic eigensolutions. In every case, non-target solutions would only be approximately recovered with an accuracy driven by the contents of the initial set of vectors. As stated above, the whole set of mono-harmonic eigensolutions may be out of reach. The reduced model however allows to recover the graph displaying the frequencies versus the dominating Fourier harmonic coefficient in the response identified from the post-processing of the generalized modeshapes. The post-processing procedures are described in the course of the considered examples.

Section 2.1 is a partial review of the literature about model reduction applied to bladed disks and rotors. Focus is set first on mono-disk reduced-order models mainly designed for mistuning analyses. The few existing reduction procedures devoted to multi-stage assemblies are then presented.

The four examples used to illustrate the methodologies are introduced in Sec. 2.2. Apart from a small but representative academic rotor model, all examples are large scale industrial models provided by Snecma.

In Sec. 2.3, one aims to present the tools included within the subsequent reduction techniques. Sector substructuring is first detailed as it is at the root of these methods. Orthogonalization procedures to turn the target vector sets into the bases of sector kinematic subspaces are then described. Finally, specific procedures devoted to the recovery of the generalized modeshapes in the case of large scale finite element meshes are considered.

Section 2.4 presents the Common Interface Component Reduction technique suited to the study of a single bladed disk with a straightforward integration of mistuning. The capabilities of this method are demonstrated with the industrial model of an integrally bladed disk.

The Disjoint Component Reduction technique is presented in Sec. 2.5. It is specially designed to deal with multi-stage bladed rotors. The main difference with the previous method is that it relies on inter-sector and inter-disk volumic interfaces to get rid of the continuity conditions that constrain the topology of the sector kinematic bases. This method is extensively tested on the proposed academic model. Its ability to handle large scale industrial models is demonstrated with the following examples.

Note finally that this chapter aims to clarify and extend material already published. Section 2.3 is quite similar to [Sternchüss and Balmès 2006], whereas Sec. 2.4 develops the orthogonalization and the reduction processes summarized in [Sternchüss and Balmès 2007a, Sternchüss and Balmès 2007b, Sternchüss et al. 2008, Sternchüss et al. 2009].

State of the art

This section aims to recall the existing reduction techniques suited to the study of single disks or multi-stage assemblies.

Purely mono-disk analysis

Many of the recently proposed techniques are only devoted to mono-disk analyses, with a clear emphasis on mistuning capabilities.

A classical modal analysis approach exists where the modes of the mistuned structure are sought as a linear combination of those of the underlying tuned structure. It is called "Mistuning Projection Method" [Bladh et al. 2001a]. [Yang and Griffin 2001] proposed a generalization of this method by seeking the approximate solution in a subspace spanned by a "Subset of Nominal Modes" (SNM), where the nominal model can be the underlying tuned configuration (it is not mandatory). The stresses in the structure are then also a linear combination of the stresses in the underlying system. The nominal modes are then found exactly with the reduced-order model. However, such methods are accurate only for a small amount of mistuning, in other words if the true kinematic subspace is not too different from the approximate one, even though it gives similar results to those returned by Component Mode Synthesis (CMS) methods in that very case [Moyroud et al. 2002]. The SNM method was then extended by [Feiner and Griffin 2002] to the case of isolated families of blade-dominated modes. In this approach called the "Fundamental Mistuning Model" (FMM), the subspace of tuned modes is supplemented by the spatial DFT of the deviations from the tuned blade frequencies. Nevertheless, the assumption of an isolated family fails in regions of high modal density. A variant of this method (FMM ID) allows to determine experimentally the origins of mistuning and to measure the natural frequencies of the tuned system [Feiner andGriffin 2004a, Feiner andGriffin 2004b].

[ Petrov et al. 2002] proposed a method that includes mistuning in a particular matrix built from (a) the transfer matrix of the tuned system and (b) a perturbation in the dynamic stiffness matrix of the system. The subsequent formulation relies on a partition of the DOF set into two families: active DOF to which mistuning is assumed applied and passive to which it is not. This method was applied to a turbine disk with non-linear contacts [START_REF] Petrov | [END_REF]].

All the techniques described above require that the amount of mistuning is small enough so that the accuracy is good even with a small set of tuned modes. But as mistuning increases, the number of tuned modes required increases dramatically until it reaches a point where the reduction becomes useless. [START_REF] Lim | [END_REF]] introduced a modal acceleration technique based on the static or quasi-static flexibilities (the latter resulting from a frequency shift of the static modes) to improve the accuracy of the forced response computations in the case of large mistuning. Such a technique is indeed well suited to deal with large changes in the geometry of one or more blades (due to FOD or blade release, for instance).

Component Mode Synthesis methods also exist to address this type of problems. They consist in partitioning the main structure into smaller components separated by interfaces. To each substructure is attached an approximate kinematic subspace spanned by a set of its modeshapes in which are found its interior motion, its interface motion, its rigid body motion... CMS methods are characterized by the boundary conditions used to compute the modes of the substructures: fixed interfaces [Craig Jr. andBampton 1968, Craig Jr 1987], free interfaces [Mac-Neal 1971, Rubin 1975, Craig and Chang 1977, Craig Jr 1987], hybrid cases [MacNeal 1971] or loaded interfaces [Benfield and Hruda 1971].

[ Yang and Griffin 1997] designed a technique that partitions a bladed disk into two substructures. The first one contains the disk, the second one contains all the blades. The kinematic subspace of the disk is built from a set of modes computed assuming that the interface between the disk and the blades is free. The kinematic subspace of the blades is built assuming that the interface between the disk and the blades has rigid body motion, it is then supplemented by the cantilevered blade modes. This method tends to overestimate the frequencies of the tuned system since the rigid body assumption overconstrains it. It also fails rapidly as the modal density increases.

[ Bladh et al. 2001a, Bladh et al. 2001b] used the same type of substructuring coupled to a Craig-Bampton reduction technique [Craig Jr. and Bampton 1968] where the coupling interface is the disk/blade interface to compute the fixed interface and constraint modes. The disk part is assumed cyclically symmetric, hence its fixed interface modes are computed with the classical one-/two-sector approach. This reduced-order model could remain quite large as is and a secondary modal analysis is performed, leading to the so-called "Secondary Modal Analysis Reduction Technique" (SMART). [Moyroud et al. 2002] employed the classical Craig-Bampton technique in a sector substructuring technique, similar to that implemented within [START_REF][END_REF]]. [Bladh et al. 2002] also used this substructuring technique with a sector subspace spanned by cantilevered blade modes and modes of the sector with a massless blade to get rid of the constraint modes. [Ratier 2004] extended this technique by reducing the sector model with the help from interface modes.

[Capiez-Lernout and Soize 2004] also employed a Craig-Bampton technique for the blade but kept the full set of interior DOF of the disk. The reduced-order model remains quite large and a second reduction is performed by considering a reduced disk model with the disk/blade interface being loaded [Capiez-Lernout et al. 2005].

[ Tran 2001] made a comparative review of many of the methods described above but in the case of a sector substructuring with inter-sector interface modes.

Multi-stage analysis

Whereas numerous methods exist to deal with single-disk analysis, only a few methods are available to deal with the multi-stage case.

As far as the author knows, the first paper to consider multi-stage coupling effects is that of [Bladh et al. 2003]. From the study of a simple but representative 3D assembly of two bladed disks, the following conclusions were drawn:

• a sort of "multi-stage mistuning" occurs because of the non matching number of blades,

• the choice of the inter-stage boundary conditions is critical if one wants to use separate reduction techniques for the individual disks and then assemble the two reduced-order models.

From these observations, the same authors built a reduced-order model of a rotor with the help from the "Component Mistuning Method" (CMM) [Song et al. 2007]. A kinematic subspace spanned by a subset of single-stage tuned modes and cantilevered blade modes is defined for each sector of each disk. The latter shapes are those to which mistuning can be applied. A secondary modal analysis allows to reduce further the number of generalized DOF so that they are attached to "tuned rotor modes" (provided that both disks are tuned). Such a method enables to

• sort the modes according to their localization, derived from the distribution of the strain energy across the rotor, even returning the dominating harmonic.

• to quantify the coupling ratio between the disk and the blades from the repartition of strain energy inside a given disk.

• to identify the mistuning ratio of each blade from a minimal set of sensors, even when that ratio is large (10%).

Similarly, [Sinha 2007] proposed a method where the kinematic subspace of the whole assembly is spanned by tuned modes of the individual disks. The reduced-order model is built from a "Modified Modal Domain Analysis" (MMDA) which allows to assemble the reduced matrices of the assembly from those of the initial sectors. This method has not been tested on a realistic finite element model yet.

Description of the numerical examples

This section aims to briefly present the numerical examples considered to illustrate the two reduction techniques described in Secs. 2.4 and 2.5.

Industrial model of a single disk [HEBD]

The first model is an integrally bladed disk provided by Snecma. It is referred to as the HEBD model in the following. The characteristics of the full finite element model are presented in Tab. 2.1. It can be divided in 23 bladed sectors whose individual mesh is described in the same table. The material is an alloy of titanium. The outer radii at the edge of the disk and at the tip of the blades are around 110 mm and 140 mm, therefore the disk is quite stiff.

Academic model of a bladed rotor [ACAD]

Preliminary analyses of multi-stage assemblies required a small but representative model of a integrally bladed rotor. To address this issue, the model displayed in Fig. 2.2, inspired by the model in [Bladh et al. 2001a], was designed. It is referred to as the ACAD model in the following sections. The radii at the roots of the blades are 200 mm and 205 mm respectively, the radii at their tips are 280 mm and 285 mm respectively. Their thickness varies from 6 mm at the root to 2 mm at the top. Their widths are 43 mm and 46 mm at the top. The characteristics reported in Tab. are those of the basic model. Indeed, its mesh is parametrable in the sense that it can be made denser or coarser at will, with either first or second order elements. The ring results from a Delaunay tessellation between the nodes of the two matching flanges of disks D 1 and D 2 . The material is titanium whose properties are set to E 120 GPa, ν 0.33 and ρ 4700 kg.m ¡3 in the material library of SDT [Balmès et al. 2007].

Industrial model of a bladed rotor [3SHP & FSHP]

Finally, the full model of an industrial rotor is considered, it is depicted in Fig. 2.3. It is composed of three integrally bladed disks (blisks), three disks (without their blades mounted), a conical drum, a sealing labyrinth and a turbine (without its blades, represented by concentrated masses). For reasons of confidentiality, the mesh cannot be displayed. tessellation, which are 4 node tetrahedrons. The size of the rings, when compared to that of the disks, is not sufficient so that their first order mesh has a significant impact on the overall solution.

Since the nominal geometries of the sectors of two connected stages are such that they are bonded together, a pre-processing of the geometries is required to slightly translate the nodes of the flanges of disks D 1 and D 3 so that the inter-disk rings can fill the empty spaces created between the flange of disk D 1 (resp. disk D 3 ) and the rim of disk D 2 , such as depicted in Fig. 2.4. Two sets of nodes have to be moved for each flange:

• the first set is translated longitudinally,

• the second set is translated radially. This operation does not affect the integrity of the finite element mesh of such parts. Notice that the bolting system is not represented in the model. The minimal set of elements to describe the geometry of the full assembly is presented in Fig. 2.5, since the bladed disks can be generated by circular repetitions of their initial sector. This allows significant memory savings because only this minimal geometry is stored. Note that without further information about the angular shifts between adjacent stages, the initial sector models of the three integrally bladed disks are circumferentially set up so that the blades are in line along e z . Albeit interesting for real assemblies, it is beyond the scope of this research to carry out studies on the influence of the relative angular positions of the initial sectors on the modal characteristics of the rotor.

Common methodologies

Substructuring into superelements

The reduction techniques described in Secs. 2.4 and 2.5 are Rayleigh-Ritz methods. They generate approximations of the kinematics of the response by combining solutions of simple problems, e.g. modes or static shapes. Linear combinations of these are assumed shapes for an approximation subspace. If the initial shapes are not independent, a basis generation procedure is needed as will be discussed in 2.3.2. The model equations are reduced by a congruent transformation that involves this basis. The number of DOF of the reduced model then corresponds to the subspace dimension.

The proposed methods are component-based in the sense that the reduction is performed for parts of the structure (in the present case, sectors of a disk). The reduced component models are similar to the pair elementary matrix-physical DOF for a finite element, hence the common denomination of superelements. The concept of superelement is illustrated in Fig. 2.6. Note that this figure exactly holds for the methodology described in Sec. 2.4 but is slightly different from the technique exposed in Sec. 2.5. Proper figures are included in these two sections. Besides, only interface DOF are displayed in Fig. 2.6b, but as will be discussed in Sec. 2.3.2, superelements may also have purely interior DOF. A classic property of Rayleigh-Ritz methods is that if exact solutions are used to generate the reduction subspace, the reduced model is exact. As detailed in Chap. 1, cyclic symmetry solutions are exact modes for tuned disks but approximate eigensolutions for rotor made of tuned disks. Using the restriction of these solutions to sector S d,s to reduce its model will thus generate a reduced model that will be:

• exact for the retained modes and approximate for others in the case of a tuned disk,

• exact for the retained mono-harmonic eigensolutions and approximate for others and for modes in the case of a rotor made of tuned disks.

Considering the graph that displays frequencies versus Fourier harmonic coefficient δ shown in Figs. 1.5 and 1.19, one retains the following shapes:

(i.) a set of mono-harmonic solutions r RepΦ d,s rcyc q ImpΦ d,s rcyc q s whose harmonics δ are selected for a particular reason (harmonic contents of the external excitations, ...). These shapes are identical for all the sectors of a given disk. They will account for either disk-dominated motion and multi-stage coupling, if applicable, or blade-dominated motion according to their position in the graph f pδq, see Fig. 1.5 and 1.19. (ii.) if needed, a set of eigenmodes rΦ d,s f ix s of this sector with its interfaces fixed. In mono-disk analyses, only the inter-sector interfaces are fixed, they are sometimes referred to as the left and right interfaces. In multi-stage analyses, the inter-disk interfaces, which are the interfaces between the disks and the volumic rings, are also fixed. The generic denomination of fixed interface modes is used throughout the rest of this chapter and no distinction is made between the two cases without loss of generality. Fixed interface modes can differ from a sector to another. They will account for blade-dominated motion.

Remark: fixed interface modes versus high δ solutions

High δ shapes would be nearly identical to fixed interface modes since they tend asymptotically towards the latter shapes. But this occurs disk by disk due to aliasing: any mono-harmonic eigensolution with a high δ relative to a given disk may fold to a lower δ relative to another disk once half its number of blades has been exceeded. This situation may lead to disk-dominated motion on one side and perturbation of the blade-dominated motion on the other side. Fixed interface modes are then preferred since they are solutions of individual disk problems. However, to keep the methodology developed for cyclic eigensolutions, such shapes can be considered as mono-harmonic eigensolutions with δ 0 with only the right interface fixed and the inter-sector interfaces also fixed, when multi-stage assemblies are considered.

The initial set of assumed vectors considered to generate the reduction will thus always take one of the two forms

rT d,s init s r RepΦ d,s rcyc q ImpΦ d,s rcyc q Φ d,s f ix s, rT d,s init s r Φ d,s f ix RepΦ d,s
rcyc q ImpΦ d,s rcyc q s.

(2.1)

A first difficulty is that such sets of vectors are often degenerate. A trivial case of dependence is found if one retains both modes of pairs found for δ ¡ 0, since it is known that the second model can be found by shifting the Fourier transform of the modeshape. Less obvious cases are found when retaining vectors from multiple δ. The independence of vectors may then be lost for the considered numerical precision. Methods to deal with this problem are discussed in 2.3.2.

A second difficulty is that coupling of components within a Rayleigh-Ritz procedure requires the assumed shapes to be continuous (kinematically admissible). A first method, described in 2.4, takes a basis of subspace generated by rT d,s init s and completes it to obtain left and right generalized DOF whose equality across sector interfaces implies continuity at all points of the unreduced model.

The first method generates bases that are possibly larger than the dimension of rT d,s init s and it is difficult to generalize to the case of multi-stage shafts where the surfaces between two disks differ from sector to sector. A second approach completely bypasses the difficulty of enforcing continuity by considering reduced models that are disjoint. Components are then coupled by reducing the model of an interface rather than enforcing inter-component continuity.

Both methods are illustrated using the various applications introduced in Sec. 2.2.

Assumed vectors and reduction bases

In the considered applications, one uses reference computations to generate sets of assumed shapes. When keeping modes of multiple target Fourier coefficients (or target rotation speeds, see Chap. 3), the sets of generating vectors kept in rT d,s init s are not bases in the sense that they are not full rank. Using non full rank reduction bases generates models that are singular or at least poorly conditioned. Procedures to generate bases from large vector sets are thus critical for the proposed reduction methods.

In linear algebra, the classical approach to generate a well-conditioned basis for a subspace is to generate an orthonormal basis of that subspace. The key issues in the process are the way to guarantee orthogonality in the presence of numerical round-off errors and the criterion used to eliminate dependent vectors.

The singular value decomposition (see [START_REF] Golub | [END_REF] for example) is the most robust approach to generate an orthonormal basis. Using a Euclidean norm on a vector of degrees of freedom is however not particularly suited for applications in mechanics. It was shown in [Balmès 1996a] that solving for the lowest frequency eigenvalues of the reduced mass and stiffness matrices gives a proper extension of the singular value decomposition to mechanics. The result of the decomposition is of the form

rTs rUsrΣsrV s T (2.2)
where rΣs is diagonal and ranks the contributions of the orthonormal vectors in rUs. A basis of the subspace is thus obtained by keeping the columns tU j u corresponding to significant singular values (such that σ j ¥ ε σ 1 ).

A proper use of the SVD requires the columns of rTs to have comparable initial scaling. This is typically achieved by initial scaling of each of the columns. In the proposed applications, the bases are composed of modes that are typically mass normalized. No rescaling is thus needed.

Gram-Schmidt methods (see Appendix 2 of [Boiteau and Nicolas 2007] for example) provide a classical procedure to generate a basis from a set of vectors. They are however very sensitive to numerical round-off. Modified Gram-Schmidt methods do a recursive orthogonalization which alleviates the numerical problem but do not introduce a mechanism to eliminate vectors. Iterative modified Gram-Schmidt (IMGS) provides a procedure that gives a robust mechanism to eliminate vectors that are not independent.

In Gram-Schmidt methods, vectors are considered dependent if the norm of the part orthogonal to the currently retained vectors is too small. Practically, if t p T j u is the part of tT j u orthogonal to tT k u with k j, the vector is truncated if

|| p T j || 2 ε ||T j || 2 .
The use of orthonormal bases in both SVD and Gram-Schmidt methods is motivated by the need to guarantee the independence of base vectors. Another approach to achieve the same result is to ensure that each vector is non zero at a different location. This is the classical approach of finite element models where shape functions are equal to one at the location associated with their DOF and zero at the locations of other DOF. This idea motivated the Independent Maximum Sequence (IMS) algorithm presented in [START_REF] Balmès | Orthogonal Maximum Sequence Sensor Placements Algorithms for modal tests, expansion and visibility[END_REF]] for sensor placement and used here to generate bases.

The IMS algorithm iteratively defines a series of independent vectors generating the target subspace as follows. At step i,

• one seeks in vector tT i u the value of the DOF m with maximum absolute amplitude. If this amplitude is too small (a typical value of 10 ¡8 ¢ maxpT 0 q), the vector is eliminated.

• the vector is normalized in the sense of norm ||.|| V t r T i u tT i u T i pmq tT i u ||T i || V ,
• further vectors tT j u with j ¡ i are made independent from tT i u by making the response at DOF m equal to zero t r T j u tT j u ¡ t r T i u T j pmq.

This step is not an orthogonalization but also guarantees independence too.

This can be illustrated as follows:

(i.) at step 0, start with the four shapes in Fig. 2.7a. The first sensor is placed where the amplitude of the first shape is maximum, that is to say at DOF 128z where the deflection of the blade in a bending motion appears to be maximum. To reach step 1, the three remaining shapes are processed such that their amplitude at this DOF is zero.

(ii.) at step 1, the second sensor is placed at DOF 57z, that is to say where the amplitude of the blade in its torsional motion is maximum. The two remaining shapes are processed such that their amplitude at this DOF is zero. (iii.) at step 2, the third sensor is placed where the deflection of the disk in its bending motion is maximum.

The algorithm iterates until the last vector is processed.

This procedure generates vectors that are obviously independent (the maximum of tT j u is reached at a position where T j 0 for j ¡ i) and normalized (the maximum of all vectors is 1).

The elimination mechanism is well suited if the initial vectors are scaled in a similar fashion, as is the case in the proposed methods where the generating vectors are mass-normalized modes. The cost of finding the maximum response is much lower than that of computing a scalar product, the procedure thus costs significantly less than Gram-Schmidt methods and does not require any iteration for round-off problems.

A significant motivation for the use of the IMS algorithm in this work is its ability to generate bases that are strictly zero on a subset of DOF. The procedure is as follows. Let the full DOF set be partitioned into two sets a and b. The modified IMS algorithm proceeds iteratively in three phases. One first reorders vectors by the maximum amplitude within DOF set a. One then selects a basis of the restriction of tTu to DOF set a. To do so, at step i,

• one seeks in vector T i the value of the DOF m with maximum absolute amplitude within DOF set a.

• if ||T i || V ε 1 ||T 0 || V ,
the vector is considered dependent and eliminated: linear combina- tions with vectors of lower index lead to a null vector.

• if T i pmq ε 2 ||T i || V , the vector is considered null within set a. tT i,a u is set strictly to t0u and the vector is set aside for the third phase.

• further vectors tT j u with j ¡ i are made independent from tT i u by making the response at DOF m equal to zero t p T j u tT j u ¡ t p T i u T j pmq.

In the last phase, vectors set aside as having a strictly null contribution on DOF set a can be made independent using the classical IMS approach.

The procedure thus generates a basis which guarantees vector independence by means of blocks of zeros. Like all methods with elimination, the results are somewhat dependent on the tolerances (values ε i ). While the issue should be considered in further detail, the results presented use ε 1 10 ¡8 and ε 2 10 ¡4 .

Partial recovery and post-processing

For very dense finite models, expanding the generalized modeshapes to the full set of DOF would be a very expensive numerical effort. In the present applications however, the motion is to be restored superelement by superelement, or equivalently sector by sector. A procedure of partial recovery of the modeshapes has been set up and is detailed further on. A subset γ of physical DOF can be chosen for each sector S d,s so that the motion is recovered to this subset only thanks to

tq d,s γ u rT d,s M,γ s tq d,s M u, (2.3)
where rT d,s M,γ s is the restriction of rT d,s M s to subset γ. This subset may vary from a sector to another. It may not necessarily derive from the true finite element mesh. Two typical choices are displayed in Fig. 2.8 and are:

• a pseudo-mesh based on nodes of the finite element model, connected by vizualisation elements. This gives only a quick insight into the motion, but allows to determine what kind of blade motion is involved, provided that the set of blade nodes is large enough to recognize it. Nevertheless, it provides a wireframe to perform test/analysis correlation.

• a subset of finite elements, which can be used for strain, stress and energy post-processing. In the case of assemblies of bladed disks, DOF set γ for partial recovery is generally chosen:

• either by a spatial Fourier analysis of the generalized displacement or the strain energy associated with each superelement. It provides useful information on where the strain energy is located (which disk and which sector in that disk). This leads to recoveries to full sector models. Further element selections in sector models are then done according to (a) the motion in the sector and (b) the user's knowledge. For instance, blade elements are selected according to its motion: the first bending motion induces high stress levels at the foot of the blade, the first torsion induces high stress levels in the vicinity of the leading and trailing edges at half the length of the blades and so on...

• or if the generalized displacement appears to be mono-harmonic, a procedure of matching mono-harmonic solutions to generalized modes (not described here) provides useful information where to look for stresses or energies.

Common interface component reduction [CICR]

The first reduction technique presented in this chapter is devoted to the study of single bladed disks. It is exactly summarized by Fig. 2.6b in the sense that each bladed sector is condensed into a single superelement connected to its neighbours by the physical inter-sector interfaces.

Basis of modes

To allow the assembly of the sector superelements into a disk model, the proposed method combines matching left and right interface motion and interior motion for which the interfaces are fixed, as illustrated in Fig. 2.9. (2.4)

The first step of the procedure is to build the subspace generated on sector 0 by the target modes (see Sec. 2.3). A basis of this subspace is built through the following procedure: (i.) One considers only the restriction of the initial set of vectors rT s init s (see Sec. 2.3 for its content) on both interfaces with the help from observation matrices rc r s and rc l s (see

Sec. 1.1.4) rT s rl s c r T s init c l T s init & .
(2.5)

If fixed interface modes are present in the initial set of vectors, they vanish naturally at this step. A SVD of rT s rl s is performed to make these vectors linearly independent (see

Sec. 2.3.2) rT s rl s rU Σ V s . (2.6)
rΣs is a diagonal matrix that contains the singular values λ j . The linear independence is guaranteed modulo the tolerance ε applied to truncate the set of singular values, leading to N λ kept vectors where

N λ maxtj, λ j {λ 1 ¡ εu.
(2.7)

The N λ kept columns of U are a basis of the image of rT s rl s, or in other words a basis of the motion restricted to interfaces I s r and I s l . Then, from basis rT s init s, one defines rT s c s

T s init v j¡N λ Φ s f ix $ , rT s i s T s init v j¤N λ $ .
(2.8) rT s c s contains both the initial vectors rejected by the SVD and fixed interface modeshapes.

An additional zero condition is imposed at the interfaces in the first equation to ensure that rT s c s will really be a basis for interior modes with both interfaces fixed, in other words a basis for modes that correspond to diagram (c) of Fig. 2.9 with rT s c s |I s r I s l r0s .

(2.9) rT s i s and rT s c s are orthonormalized with respect to the stiffness rK s s of sector S s through a standard Gramm-Schmidt procedure

rT s i s rT s i s ¡ rT s c t K s T s c s ¡1 rT s c s rT s c t K s T s i s.
(2.10)

This leads to bases whose general forms are

rT s c s ! 0 T s c,c 0 ( ) , rT s i s ! T s i,r T s i,c T s i,l (
) .

(2.11) (ii.) the following step is to build the bases of the motion at the interfaces from rT i s. This set of vectors is supplemented such that

rT s i s T s i1 T s i2 $ , (2.12) with rT s i1 s ! T s i,r T s i,c T s i,l ( ) , rT s i2 s ! T s i,r ¡K s c,c ¡1 K s c,r T s i,r 0 ( ) .
(2.13) rT s i2 s contains the right interface vectors, their static recovery to the complementary domain and no motion at the left interface.

(iii.) Right and left interface vectors are matched modulo θ, which is the rotation of angle α 2 π{N s . Introducing the following matrix rB s s

0 I ¡ c r r T i2 © ¡1 ¡ c l r T i1 © ¡I ' , ( 2
$ ! θpT s i,l q ¡K s c,c ¡1 K s c,r θpT s i,l q 0 ( ) , T s S,l $ ! 0 T s i,c K s c,c ¡1 K s c,r θpT s i,l q T s i,l ( ) .
(2.16) (iv.) With rT s S,c s rT s c s, the final basis is rT s S s

T s S,r T s S,c T s S,l $ .
(2.17)

This basis has the particular shape described by Fig. 2.9. rT s S,r s is a basis of modes with the right interface free and the left interface fixed. rT s S,c s is a basis of interior modes with both interfaces fixed. rT s S,l s is a basis of modes with the left interface free and the right interface fixed.

Note that when the disk is tuned rT s S s rT 0 S s, ds v0, N s ¡ 1w. Otherwise, rT s S,r s rT 0 S,r s and rT s S,l s rT 0 S,l s, ds v0, N s ¡ 1w, whatever the degree of mistuning of each sector S s . The latter drives the contents of rT s S,c s which may vary from sector to sector. When it happens that mistuning is small enough so that the blade shapes remain unchanged, the assumption rT s S,c s rT 0 S,c s, ds v0, N s ¡ 1w remains valid, mistuning then acts on the reduced stiffness matrix.

Reduced matrices

Any finite element matrix rA s s relative to sector S s is projected onto the approximate subspace so that rA s S s rT s S t A s T s S s.

(2.18) rA s S s takes a particular form depending on whether rA s s is rM s s or rK s s since rT s S,c

t M s T s S,c s rIs , rT s S,c t K s T s S,r s T t Sc K s,0 T Sl $ r0s , rT s S,c t K s T s S,c s Λ s S,c $ z ¡ ω s cj © 2 z & . (2.19)
where ω s cj is the j th natural frequency of sector S s with fixed interfaces. As a consequence, rM s s and rK s s have then the remarkable forms

rM s S s " ! T s S,r t M s T s S,r T s S,r t M s T s S,c T s S,r t M s T s S,l T s S,c t M s T s S,r I T s S,c t M s T s S,l T s S,l t M s T s S,r T s S,l t M s T s S,c T s S,l t M s T s S,l ( 0 )
,

rK s S s ! T s S,r t K s T s S,r 0 T s S,r t K s T s S,l 0 Λ s S,c 0 T s S,l t K s T s S,r 0 T s S,l t K s T s S,l (
) .

(

2.20)

There is no stiffness coupling between the interface and the interior generalized DOF. The coupling is fully integrated in the non-zero terms of the reduced mass matrix. The matrices displayed in Fig. 2.10 are obtained with 2 ¢ 16 interface vectors and 9 interior vectors. Right and left interface vectors are built so that they are matched. This ensures that the inter-sector continuity condition (1.14) is equivalent to stating that left generalized DOF of one sector are equal to right generalized DOF of the next one rc l T s S,l s rc r T s 1 S,r s ðñ tq s S,l u tq s 1 S,r u, (2.21) where, as stated in Sec. 2.4.2, rT s S,l s rT 0 S,l s and rT s S,r s rT 0 S,r s, ds v0, N s ¡ 1w, even in the case of a mistuned disk. This condition makes the assembly of a disk model with respect to the generalized coordinates very straightforward, leading to block diagonal matrices obtained by the assembly of those of the sector superelements, as depicted in Fig. 2.11. It also provides a very convenient way to have a quick insight into the spatial distribution of the motion by defining a pseudo-mesh that represents the generalized degrees of freedom. Such a representation is depicted in Fig. 2.12. where r∆ s s is a diagonal matrix containing the mistuning factors. For a given sector, these factors can be different from a mode to another.

Verification in the case of a tuned disk

Before any consideration of mistuning, it was necessary to verify the reduction technique on a tuned disk, in particular to be sure that the modes used to build the reduction basis are obtained exactly with the reduced model. A reduced model of the HEBD is built using a basis composed of 2 ¢ 16 interface modes built from 3 single or paired modes with δ v0, 3w and a frequency in the band r0, 5s, supplemented by 9 modes of the sector with fixed interfaces.

Table 2.5 summarizes the number of DOF of the different models of this bladed disk. Figure 2.13 illustrates the contents of the reduction basis. The modes obtained with the reduced-order model are compared to the modes given by a classical cyclic symmetry approach. The resulting graph f vs. δ and the corresponding relative error in frequency for the single or paired modes with δ v0, 11w in the frequency band r0, 15s are depicted in Fig. 2.14. As expected, this error is zero for the targeted values of δ and stays under 0.7% for the non-target modes in a frequency band three times wider than that of the target modes. One notices that the use of sector modes with fixed interfaces overstimates by around 0.1% the frequencies of cyclic modes with high values of δ. The correlation between the reduced modes and the reference modes is computed by means of the MACM criterion introduced in Sec. 1.1.7 and is very good: beyond 99.99% for the target modes and no less than 93% for other modes in the frequency band r0, 15s. Maximum accuracy is obtained for modes with either small δ at low frequencies or high δ as they are very close to that of the reduction basis. As a result, modeshapes obtained with this tuned reduced model are very close to the exact tuned modes. Typical modeshapes are shown in Fig. 2.15 before and after their recovery to the full disk. They clearly exhibit similar patterns under these two forms, which demonstrates the convenience of the reduced representation. Recall that the reduced-order model is periodic, therefore all superelements are identical. The spatial harmonic content of the generalized modeshapes can be analyzed either from the reduced or recovered modeshapes. The Fourier spectra of the already considered solutions are reported in Fig. 2.16 and exhibit a single harmonic, as expected.

Application to mistuning

Mistuning is taken into account as a distribution of frequencies of fixed interface modes found in the diagonal of rK S s. A sample distribution is shown in Fig. 2.17. The reduced stiffness matrix is no longer made of identical blocks along its diagonal. As a result, the generalized modes are spatially multi-harmonic. This is illustrated by Figs. indicates that this mode only involves the third interior DOF and the visual inspection of the recovered modeshape confirms that it corresponds to the first torsion of the blade. However, the spatial spectrum of this modeshape is not strictly a Dirac's comb and the energy is not totally confined to this blade. Indeed, a few blades exhibit a small but non-zero generalized/recovered amplitude. Conversely, Fig. 2.19 presents a degenerate tuned mode with a contribution of all blades with a quite distinguishable pattern confirmed by the spatial spectrum.

The response of the disk to a stationary excitation with δ 3 is obtained directly from the reduced-order model. Figure 2.20 shows the maximum amplitude of the displacement recovered ) 2, (c) 3, (d) 5, (e) 6, 7, 8 and(f) 9 at the blade edges (tip, leading and trailing edges). The modal damping ratio is low enough so that peak scattering can be seen. The amplitude of the mistuned response is not necessarily higher than its tuned counterpart even though increases in the amplitude levels are observed.

To conclude, the aim of this mistuning analysis is not to characterize this specific disk but to demonstrate that the method is suited to this type of analyses. The direct integration of frequency mistuning within the reduced stiffness provides a very useful means to perform statistical analyses on a population of bladed disks. For each sample distribution of frequencies, careful post-processing is required to obtain the stress levels within the blades for HCF analyses, especially when a highly localized response is observed.

Disjoint component reduction [DJCR]

The main problem of the CICR technique is that two adjacent superelements must share exactly the same set of interface DOF, hence it requires numerous vector operations that can be very time and memory consuming for large sets of physical DOF. To avoid this problem, a second reduction technique is introduced where each bladed sector is divided into two superelements:

(i.) the first one is associated with all finite elements except those with nodes belonging to the left interface; this type of superelements are called (intra-)sector superelements in the rest of this section, (ii.) the second one contains the elements that have been removed in (i.); this type of superelements are referred to as inter-sector superelements.

Approximate subspace in a mono-disk approach

The proposed methodology is first introduced within a mono-disk approach. As stated at the beginning of this section, each bladed sector is condensed into a sector superelement (white region with blue edges in Fig. 2.21) that is connected to other sectors thanks to inter-sector slices of finite elements (green region with black border in Fig. 2.21). Efforts are concentrated upon the search for an appropriate subspace to approximate motion within the sectors superelements. The interface (green inter-sector portion) is later reduced by considering the trace of the motion of the adjacent sectors to the right and left interfaces of the inter-sector slices. (2.23)

The reduction process aims to build an approximate kinematic subspace of sector S s in order to

recover tq s u by tq s u rT s M s tq s M u ! T s M,r T s M,c T s M,l ( ) tq s M u, (2.24)
where tq s M u contains the generalized degrees of freedom and rT s M s is a basis of the subspace generated by rT s init s and partitioned according to Eqn. (2.23). Any finite element matrix rA s s relative to sector S s is then projected onto the subspace spanned by rT s s thanks to the relation

rA s M s rT s M t A s T s M s.
(2.25)

According to Eqn. (2.23), the slice C s that connects sectors S s and S s 1 couple tq s l u and tq s 1 r u, so that any matrix rA s C s relative to this slice has to be reduced on the subspace generated by rT s M,l s and rT s 1 M,r s such that

rA s M,C s T s M,l 0 0 T s 1 M,r & t rA s C s T s M,l 0 0 T s 1 M,r & . ( 2

.26)

Columns in rT s 1 M,r s are assumed to be vectors that span the subspace of motion on left interface I s 1 l (the s 1 will not be omitted for clarity). For a large set of generating vectors in rT s init s, this rectangular matrix is not necessarily full rank so that a generating basis must be built first.

SVD or Gram-Schmidt methods (see section 2.3.2) only work on rT s M,l s. The only way to ensure that the reduced basis generates the same subspace as rT s init s is to decompose the basis in two parts and apply an SVD or Gram-Schmidt on both

tq s u " ! r T s M,r 0 0 r T s M,c 0 r T s M,l ( 0 ) .
(2.27)

With the modified IMS algorithm detailed in section 2.3.2, the topology of the matrix is

tq s u " ! r T s M,r 0 T s M,c T s M,l ( 0 )
4

q s M,r q s M,c B .
(2.28) with a number of vectors that is smaller than or equal to the size of rT s init s. Even though it is never explicitly built, the subspace in which the approximate solutions of the disk problem are sought is given by

rT M s " ! T 0 M ¤ ¤ ¤ 0 . . . . . . . . . 0 ¤ ¤ ¤ T N 1 s ¡1 M ( 0
) .

(2.29)

Approximate subspace in a multi-disk approach

When it comes to multi-stage assemblies, sectors of a given disk are not only connected to each other, but with one or more sectors that belong to the adjacent disks. In this case again, each bladed sector is condensed into a sector superelement (white region with blue border in Fig. 2.23) that is connected to its neighbours thanks to inter-sector slices of finite elements (green region with black border in Fig. 2.23) and inter-disk rings R (red region).

The DOF set of sector S d,s is partitioned according to Fig. 2.23 so that tq d,s u 6 9 9 9 9 9 8 9 9 9 9 9 7

q d,s rR q d,s r q d,s R q d,s c q d,s l D G G G G G F G G G G G E .
(2.30)

The approximate kinematic subspace of sector S d,s is built so that tq d,s u is recovered by tq d,s u rT (2.33)

T d,s M,r T d,s M,R T d,s M,c T d,s M,l ( 0 0 0 0 0 0 ) . ( 2
According to Eqn. (2.30), slice C d,s that connects sectors S d,s and S d,s 1 couple tq d,s l u to tq d,s 1 rR u and tq d,s 1 r u, so that any matrix rA d,s

C s relative to this slice has to be projected on the subspace generated by rT d,s M,l s, rT d,s 1 M,rR s and rT d,s 1 M,r s such that

rA d,s M,C s " ! T d,s M,l 0 0 T d,s 1 M,rR 0 T d,s 1 M,r ( 0 )
t rA d,s C s " ! T d,s M,l 0 0 T d,s 1 M,rR 0 T d,s 1 M,r ( 0 )
rT M,R s " " " " " " " " " " " " " " " " " " " " ! T 1,0 M,rR ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 T 1,0 M,R ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ T 1,N 1 s ¡1 M,rR 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ T 1,N 1 s ¡1 M,R 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 T 2,0 M,rR ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 T 2,0 M,R ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ T 2,N 2 s ¡1 M,rR 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ T 2,N 2 s ¡1 M,R ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ) . ( 2

.35)

Any matrix rA R s relative to ring R is projected on the full subspace spanned by the vectors in rT M,R s rA R M,R s rT t M,R A R T M,R s.

(2.36)

A careful sorting of the DOF sets and a smart use of the various orthogonalization strategies described in Sec. 2.3.2 enable to limit connectivity and guarantee vector independence. One targets the typical basis topology Once more, the full subspace rT M s is not explicitly built as .38) but the reduced matrices are assembled from the superelementary matrices in two steps:

rT M s " " " " " " " " " ! T 1,0 M ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ T 1,N 1 s ¡1 M 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 T 2,0 M ¤ ¤ ¤ 0 . . . . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 0 ¤ ¤ ¤ T 2,N 2 s ¡1 M ( 0 0 0 0 0 0 0 0 0 ) , ( 2
Step 1: the physical finite element matrices of the sectors and inter-sector superelements are built sector by sector, disk by disk, thanks to Eqns. (2.33) and (2.34), leading to matrices rA d,s M s and rA d,s M,C s associated with sector S d,s and slice C d,s respectively. These matrices are assembled into reduced disk matrices, exactly like in the mono-disk approach.

Step 2: the matrices of the ring are projected onto the kinematic subspaces of its neighbouring disks with respect to the inter-disk generalized DOF thanks to Eqn. (2.36). In order to avoid handling the full matrix of rT M s, rA R M s is built iteratively block by block, or equivalently sector by sector.

Table 2.6 summarizes the memory required to store both the initial matrix and its factors with or without the explicit isolation of the DOF common to the intra-sector and right inter-sector interfaces. At this point of the development, these figures do not seem in favour of the isolation procedure whereas the matrices look sparser (see Fig. 2.25 for a model with 3 and 4 sectors). Special efforts will be paid to that point in the upcoming work to demonstrate that such a feature is a true improvement.

Application to the academic rotor [ACAD]

This reduction technique is first applied to the academic rotor model presented in Sec. 2.2.2 to study its convergence.

The reduction process leads to relatively sparse matrices with the remarkable form displayed in Fig. 2.26 for this assembly of two tuned disks with 12 and 15 blades. The diagonal blocks and four off-diagonal blocks derive from the projection of the sector matrices. The latter blocks can be divided into two groups, the first with 12 blocks and the second with 15 blocks. The inter-disk coupling is contained both in the diagonal blocks that derive from the matrices of the inter-sector joints and in the off-diagonal blocks that derive from the matrices of the ring. Apart from these remarks, the reduced matrices have no other properties which could impact their eigenvalues or eigenvectors.

Different generating sets of mono-harmonic solutions with δ v0, 5w are considered in various target frequency bands displayed in Fig. 2.27a. 250 generalized modes with a frequency f r0, 15s

Application to the 3-stage industrial model [3SHP]

A sub-model made of the 3 integrally bladed disks of the industrial rotor depicted in Fig. 2.3 is now considered. The minimal geometry composed of one sector per disk plus the rings derived from that in Fig. 2.5 is perfectly suited to the computations of the mono-harmonic eigensolutions. A classical graph of normalized frequencies versus Fourier harmonic coefficients is plotted in Fig. In particular, the recovery of the strain energy in Fig. 2.29b provides useful information since its maximum is located at the rim between disks D 1 and D 2 and not at the tip of the blades as would be expected because of the 2S1 shape. Particular attention must be paid to the fact that for a tuned disk the stress and energy distributions are similar on every sector except that levels follow a sine/cosine pattern driven by δ. Hence, the stress (resp. energy) level retrieved on sector S d,0 must be considered with care since it may not be the maximum stress (resp. energy) withstood by disk D d . This emphasizes the need for 360°disk/rotor models. A reduced model of this assembly is built from an initial set of target mono-harmonic eigensolutions with δ t0, 1, 2u and f r0, 16.13s plus fixed interface modes in the same frequency band.

This reduced model has 44 342 generalized DOF. The topology of the reduced stiffness matrix is displayed in Fig. 2.33. 2 Gb are required to store either the mass or the stiffness matrix. Once factored, any of these matrices occupies 16 Gb on the hard disk drive.

900 modes are computed to cover the entire normalized frequency band r0, 19.5s. Before restoring the modeshapes to the whole structure or a subset of physical DOF, a convenient way to describe the generalized modes is to perform a Fourier spatial analysis of the reduced modeshapes. This analysis is conducted disk by disk since their reduced models are rotationally periodic, which is obviously true for tuned disks but also for disks with a mistuning small enough so that it does not affect the sector subspaces. Mistuning is not considered here and this point is discussed further in Sec. (2.39) and is plotted in parallel to these spectra. It is mandatory to examine both graphs together to avoid hasty conclusions in terms of multi-harmonicity of modes. Indeed, a mode that looks multi-harmonic on a given disk can be so if its participation relative to that disk is very small: the analysis of an almost zero signal leads to a spectrum that shows multiple harmonics, but it does not indicate that the generalized displacement is almost zero on this disk! A physical indicator, such as the kinetic energy, should be used to avoid normalization problems between the kinematic subspaces defined on each disk. It is introduced in the following example. For instance, modes 118 ¡ 180 are clearly localized to disk D 3 according to Fig. 2.34 but seem to be multi-harmonic on disks D 1 and D 2 . The only multi-harmonic modes of the reduced model are observed on disk D 3 for δ ¡ 22 (if δ is considered as the dominating harmonic in these particular cases), as described in Fig. 2.37 where the shape of mode number 428 at f 10.23 is restored to the whole assembly and its weighted spectrum is computed for the generalized vectors by

dd t1, 2, 3u, R d tq d M u t rK d M s tq d M u tq M u t rK M s tq M u ,
max dt1,2,3u ¢ max δv0,N d s {2w § § § § § §p q d δ § § § § § § 1 1,
(2.40)

This shape seems localized to the blades, thus this shape lies in an accumulation band. Accumulations are generally not affected by multi-stage coupling and the scattering of the spectrum is quite surprising. Conversely, from the latter figures, it can be concluded that almost all modes are mono-harmonic, like that of Fig. 2.38, even for highly coupled modes, which are those most likely affected by multi-stage coupling. A possible explanation is that the asymmetry of the intermediate ring is not sufficient to strongly affect the overall motion of the bladed assembly. Full recoveries are in fact not necessary to draw the same conclusions. Actually, the partial recovery described in Sec. 2.3.3 is particularly useful here since storing the 900 generalized modes represents about 305Mb, whereas storing the 900 modes once restored would represent 53Gb. To address this issue, two different selections are considered here for post-processing:

• to display the motion, a coarse mesh connecting a few nodes in each sector is depicted in Fig. 2.39a, it includes a sufficient number of blade nodes to discriminate between the various types in their motion,

• to compute stresses and energies associated with a certain type of motion, a selection of blade and disk elements is performed for disk D 1 and is shown in Fig. 2.39b. However, a partial recovery to a viewing mesh as is done in Fig. 2.40 only provides an indication where to look for high densities of strain energy and high stress levels but lacks the possibility to indeed recover such quantities since there is no underlying mechanical content in this mesh. To do so, the selection of true finite elements displayed in Fig. 2.39b is invoked and the strain energy and the Von Mises stress distributions are recovered to that subset of elements for the generalized mode at f 15.89, which is depicted in Fig. 2.41. This figure shows that, for this very mode and this very disk, (a) the Von Mises stress is maximum on the disk flange and (b) as expected it has a high value in the middle of the blade tip due to the 2S1 shape. Without proper considerations of multi-stage coupling, it is highly possible that conclusion (a) would be missed. Besides, the mono-harmonic solution also shows a strong participation of disk D 2 and the subset of kept finite elements could be extended to that disk.

Application to the full industrial model [FSHP]

In this section, the full rotor displayed in Fig. 2.3 is considered. Mono-harmonic eigensolutions are computed in the frequency band f r0, 24.17s, if the frequencies are normalized with respect to the lowest frequency among these solutions. Because of the high amount of CPU time requested to compute them within MDNastran (see Tab. 2.8 for details), it is mandatory to compute only the target mono-harmonic eigensolutions required in the subsequent reduction process. It emphasizes the need for an a priori guess of the dynamics, which relies on the user's experience. The graph of normalized frequencies versus Fourier harmonic is plotted in Fig. 2.42. Even with this small set of solutions, similar conclusions of those drawn from Fig. 2.28 can be made.

In particular, accumulation frequencies (frequencies of single sectors with fixed interfaces) can already be seen for disk D 1 (2F, 1T), disk D 2 (1F) and disk D 3 (1F). Once more, the direct display of the complex sector shapes seems sufficient to conclude about the localization of the motion at the rotor level. The localization at the disk or sector level can be achieved by means of a partial recovery of the mono-harmonic shapes to the pseudo-mesh depicted in Fig. 2. 51 and not yet used at this step of the analysis. In particular, the global shaft modes are retrieved as well as highly localized shapes.

A subset of mono-harmonic eigensolutions with δ v0, 4w and f r0, 24.6s is used to build a reduced model. The initial set of complex vectors is made of 594 vectors distributed as reported in Tab. 2.9. Building the approximate kinematic subspace of each superelement requires 4 825 s and leads to bases whose main characteristics are gathered in Tab. 2.10. Note that the amount of memory required to store each single basis is quite large. The fact that each bladed disk is tuned (or lightly mistuned, see Sec. 2.6.2) helps to limit the total amount of memory since the recovery matrices are actually stored once and reused for all sectors of a given disk. As described in the previous example, some post-processing operations can be performed before any restitution to any set of DOF. First, the frequencies resulting from the resolution of the reduced matrix problem can be plotted so that accumulations are seen. Figure 2.48 then shows that, as expected, the presence of the intra-/inter-sector fixed interface modes force the frequencies to accumulate towards these shapes. Such a graph is a useful tool in the sense that it allowed to update the orthogonalization procedures presented in Sec. 2.3.2 to force the fixed interface modes not to be eliminated, therefore to force the accumulation of the generalized modes towards these modes. Second, the spatial harmonic content of the generalized modeshapes can be investigated for the same reason as in Sec. 2.5.4: because the reduced models of the disks are geometrically periodic. The Fourier spectrum is here determined using an energy criterion. The fraction of strain energy per harmonic per disk is considered. Indeed, by virtue of Eqn. (1.46), the norm in strain energy of tqu (physical or reduced) restricted to disk D d is the sum of the norms in strain energy of each spatial harmonic tp q δ u (more precisely of its real and imaginary parts when complex). The fraction x W d δ of the total strain energy W associated with each harmonic tp q δ u is then defined by

x W d 0 § § § § p q d 0 § § § § K d,0 ||q d || K , x W d δ § § § § Re p q d δ ¨ § § § § K d,0 § § § § Im p q d δ ¨ § § § § K d,0 ||q|| K , x W d N d s 2 § § § § § § § § p q d N d s 2 § § § § § § § § K d,0 ||q|| K . ( 2

.41)

These quantities indicate which disks contain energy and which harmonics contribute. It supplements the information about the spatial mono-or multi-harmonicity of the displacement. It also avoids the plots of two separate graphs as done in Figs. 2.34 to 2.36 . Two ways exist to represent these fractions of strain energy per harmonic:

• harmonics are grouped by disk, as done in Fig. 2.49: it clearly shows single disk motion in some frequency bands with harmonics increasing towards blade modes (accumulation), • disks are grouped by harmonic (with the appropriate number of zeros for non-present harmonics), as done in Fig. 2.50: it indicates if all disks share the same subset of spatial harmonics.

Indeed, in the normalized frequency band r2.90, 6.44s Fig. 2.49 shows that an isolated accumulation is observed on disk D 3 and 2.50 shows that all coupled modes in this band are spatially mono-harmonic with low identical values of δ on each disk (typically 0 to 2). Such observations are typical of the whole frequency band r0, 29s except that a few multi-harmonic coupled modes are observed between f 14 and f 17 (a typical shape is depicted in Fig. 2.57). Such obser- vations are similar to that in Sec. 2.5.4, probably for the same reason: the inter-disk rings are small enough so that they do not affect much multi-stage motion by scattering the spatial spectra.

The harmonic decomposition of the modeshapes is supplemented by a quick view of them partially recovered to the pseudo-mesh depicted in Fig. 2.51 built with the same constraints as in the previous example. In parallel, the grouping of harmonics disk by disk is preferably used, as it is done in Figs. 2.52 to 2.57, where the localization information derived from the bar graphs is confronted with the partial recoveries of the modeshapes to the pseudo-mesh. In all cases, the partial recovery provides the supplementary information about where the motion (or equivalently the energy) is really localized within the individual disks. The study of the fraction of strain energy per harmonic per disk of all the generalized modes in the frequency band r0, 29s demonstrates that, except for a few of them, they are spatially mono-harmonic. This indicator enables to fill the incomplete graph in Fig. 2.42 by identifying the harmonic with the highest fraction of strain energy (all disks taken together) for each mode.

The result is depicted in Fig. 2.58. Accumulations are well reproduced and, more important, aliases are not present. The identification of the accumulations is made easier with the help from Fig. 2.59 where the most participating disk is indicated for each generalized mode. The information about the most energetic disk in the region of disk-dominated motion should be considered carefully since for coupled modes, the relative size between the disks makes disk D 4 the most energetic whereas the amplitude of the displacement may be more important on the others.

This reduced model is used to compute transfer functions of the rotor when two shakers are positioned at 90°of each other, along the inward-pointing normal (¡e x and ¡e y respectively) at the lower extremity of the bearing of disk D 2 . In the following, only the first excitation is considered since similar results are obtained for the second one. Each shaker excitation concerns a single DOF of disk D 2 , hence their spatial spectrum relative to this disk is a Dirac's comb as displayed in Fig. The spectra of the generalized modes that lie in this frequency band are already plotted in Figs. 2.49 and 2.50. The contribution of the individual modes to the global response can be estimated by computing the ratio of effective stiffness associated with each mode. The effective stiffness is such that N j1 tbu t tΦ j u tΦ j u t tbu ω 2 j tbu t rK ¡1 s tbu 1, (2.42) where tbu is the controllability vector, here limited to a single DOF. The effective stiffness is the sum of the contributions κ j of each mode pω j , tΦ j uq κ j tbu t tΦ j u tΦ j u t tbu ω 2 j tbu t rK ¡1 s tbu tbu t tΦ j u tΦ j u t tbu

ω 2 j 1 °N k1 tbu J tΦ k u tΦ k u J tbu ω 2 k tbu t tΦ j u tΦ j u t tbu ω 2 j 1 °Nv k1 tbu J tΦ k u tΦ k u J tbu ω 2 k , (2.43)
since the mode set is truncated to keep only N v modes within a given frequency band (not necessarily restricted to the band in which the response is computed). Effective stiffnesses are plotted in Fig. 2.62. Generalized modes are assumed controllable when their effective stiffness exceed a given threshold (here ε 10 ¡10 ), their natural frequencies are later plotted as vertical dotted lines in all graphs. The vertical green dashed line stands for the accumulation frequency of disk D 3 that lies in this frequency band. Near-pairs of modes can clearly be seen in this figure. The transfer function between the first shaker and the blades in the pseudo-mesh of Fig. 2.51 is computed in 34 587s. Once more, a procedure of partial recovery is employed to restore the response expressed on the generalized DOF set to this blade DOF set. Figure 2.63 shows the maximum response on this DOF set stage by stage. The multivariate modal indicator function (MMIF) is added to this plot to confirm the detection of phase resonances (as it indicates when the transfer functions become purely imaginary [Géradin and Rixen 1993]).

One is also interested in the spatial spectrum of the response. Two ways are possible:

• after the recovery (full or partial), a Fourier analysis of the response at a given set of sensors is performed. This method is adapted to experimental data, that is to say when only physical measurements are available. In this study, a single physical DOF is chosen for each blade of each disk so that all resonant modes in the desired frequency bandwidth are observable at this physical DOF. This choice is made here by considering the z direction at nodes at the corner of each blade with the z coordinate being minimum, as reported in Fig. 2.64.

Besides, since interest is focused on the harmonic content at the resonances, the transfer functions are processed only at these frequencies, indicated by vertical dotted lines in Fig. 2.63. The result of the processing is displayed in Fig. 2.65. The relative amplitude of each harmonic with respect to the dominating one at each resonance is plotted. It seems that only δ 1 is involved in the response, except near f 5 where a few more harmonics contribute to the response, especially in disk D 2 .

These results are better understood if one plots in parallel the measured signal at the sensors for each resonance. This is done in Fig. 2.66 for the third and seventh peaks or paired peaks, since modes come mostly in pseudo-pairs. However, in the present case, the double peaks are not distinguishable here due to a too high damping ratio relatively to the difference between the frequencies in the near-pairs. The amplitude of the response indeed present two lobes characteristic of δ 1, when plotted against the angular location of the retained sensors.

• however, the recovery is still not mandatory to conclude since one can compute the fraction of strain energy per harmonic per disk at a given frequency. Such results are presented in Figs. 2.67b and 2.68b for peaks number 3 and 7 identified in Fig. 2.66. However, if the user wants some more insight in the operational deflection shapes, a partial recovery to the vizualisation mesh is of great help to describe the longitudinal and circumferential location of the motion, as depicted in Figs. 2.67a and 2.68a. In the present case, the external load tends to direct the motion along the direction of excitation, an information that is not directly available from the study of Figs. 2.67b and 2.68b.

Conclusions

In this chapter, two reductions techniques are proposed to deal with assemblies of bladed disks. This section aims to summarize the main features of these techniques, with a particular emphasis on the methodological aspects.

Summary of the methodology

Indeed, beyond the results obtained for various configurations (single-disks or multi-stage assemblies), it is mandatory to insist on the methodology employed to deal with such problems. The two methods presented in this chapter are used in a very similar way described in Fig. 2.69. The only difference lies in the input sets of target vectors.

Future work

The two reduction techniques introduced in this chapter pave the way for efficient analyses of multi-stage bladed rotors. However, a few further developments deserve to be considered:

• the orthogonalization strategy used in the DJCR technique must be updated so that there is as little coupling as possible in the reduced matrices. Since the bandwidth drives the size of the factorized matrix, it is mandatory to keep it as small as possible, that is to say to have as few off-diagonal terms of inter-disk coupling as possible. Moreover, since the inter-disk rings are small compared to the disks, it could be interesting to investigate how neglecting their masses could impact the resolution of the reduced matrix problem.

Besides, such a reduction technique should be included into a subspace iteration strategy so that the initial set of vectors is supplemented in order to minimize an identified residue.

• parametric studies must be made available, including thermal and rotation effects. Thermal effects are critical in full rotor models since there exists a gradient of temperature, globally from the turbine to the front bladed disk, and locally from the blade tip to the inner radius of the rim. Rotation effects are discussed in Chap. 3.

• mistuning does not appear in the flowchart in Fig. 2.69. Both methods can be seen as a jigsaw puzzle where in absence of mistuning all pieces are identical. To introduce mistuning, one can modify the characteristics of each superelement independently, its mechanical properties (small mistuning) as well as the geometry of the underlying sector model (small and large mistuning).

When mistuning is small enough, the blade modeshapes are not significantly affected whereas their natural frequencies are. This is typically the case when one is interested to quantify the effects of machining tolerances to the dynamic behaviour of a family of disks or rotors manufactured on the same assembly line. In that very case, since the previous assumption concerns the fixed interface modes, it is possible to keep the same sector kinematic subspaces. One however slightly changes the diagonal terms corresponding to these shapes in the reduced stiffness, as it is done explicitly in Sec. 2.4.4. This is particularly convenient since the generalized DOF are kept unchanged, as displayed in Fig. 2.70a.

When the latter assumption fails for larger changes, in particular when a blade encounters an impact (during a foreign object damage event) and is strongly deformed, it is mandatory to update its modeshapes. As a result, the interior generalized DOF are changed and one must pay particular attention to keep the same set of generalized inter-sector DOF in the subspace construction phase, as depicted in Fig. 2.70b. • finally, it would be particularly interesting to confront such a numerical model to experimental data obtained on a multi-stage sample rotor so that mono-or multi-harmonic modeshapes could be observed.

3

Structures with Variable Rotation Speeds The rotation of any structure induces inertial loads that cannot be neglected when the rotation speed is high. This is all the more critical for slender structures such as blades in turboengines, especially for the fan and booster stages where long wide-chord blades are found.

Inertia leads to quasi-static deformations in such blades, mostly an extension with a subsequent detwisting. As a result, the aerodynamic performances of the blades are locally altered which affects the flow in their vicinity and consequently modifies the characteristics of the fluid-structure interaction. The decrease of the gap at the tips of the blades may also reduce the aerodynamic efficiency of the whole stage at best or lead to a contact with the casing at worst [START_REF] Legrand | [END_REF].

Indeed, rotation-induced forces are follower forces, thus the two-step procedure described in Sec. 3.1 has to be employed. Both problems rely on the same weak form which is recalled in parallel to the associated discrete form so that supplementary operators or matrices are seen. The first step computes the deformed state of the structure when it is brought from rest to a quasi-static equilibrium induced by the inertial load. This operation involves either linear or non-linear statics. The second step computes the free or forced response of the structure once this prestressed state has been reached. The effects of the supplementary terms in the matrix problem on the dynamics are then discussed. Among them, it is observed that eigenfrequencies evolve with respect to the rotation speed due to stiffening/softening effects. In addition, gyroscopic coupling also affects them as it tends to separate double frequencies (at rest) into two distinct ones which correspond to a forward-and a backward-travelling wave respectively [Jacquet-Richardet et al. 1996,Chatelet et al. 2005[START_REF] Santos | Contribution to experimental validation of linear and non-linear dynamic models for representing rotor-blade parametric coupled vibrations[END_REF]]. The evolution of frequencies is usually depicted in a graph ω j f pΩq called a "Campbell diagram" presented in details in Sec. 3.2.5. Such graphs are of particular importance to deal with engine order excitations whose frequencies depend linearly on the rotation speed of the engine (see Remark page 22). However, such supplementary matrices have to be assembled at each rotation speed of interest. It is demonstrated in Sec. 3.2 that the supplementary stiffness matrices can be replaced by a matrix quadratic polynomial of the squared rotation speed Ω. This polynomial is actually exact in certain limit cases described in this section, only an approximation otherwise. This enables to cover wide ranges of rotation speed with only three constant matrices, derived from exact stiffness matrices computed at three rotation speeds.

Even with such a convenient tool, the numerical cost to deal with various operating points in a given speed range may become prohibitive. As a result, the combination of (i) the approximation of the supplementary stiffness matrices by a quadratic polynomial of the squared rotation speed introduced in Sec. 3.2 with (ii) a multi-model approach where an approximate subspace is built from a set of exact vectors computed at three rotation speeds described in Sec. 3.3 provides a means to significantly reduce that cost. These methods are extensively tested on various models shown in that section.

Note finally that Sec. 3.1 highly supplements [Sternchüss and Balmès 2006] in the sens that the emphasis is put on the underlying finite element formulations from which the polynomial approximation is derived.

Theoretical background

As described in the introduction to this chapter and summarized in Fig. 3.1, computing the response of a given structure while in rotation requires to solve two successive problems: (i.) first, the shape associated with the static-stress state must be determined with the help from the large displacement theory, (ii.) second, for small amplitudes of vibrations, they allow to linearize the dynamic problem around the prestressed state.

In other words, what is desired in this section is to formulate the evolution problem that turns the initial configuration V p0q (the structure at rest whose geometry is assumed well known and meshed) into a first quasi-static configuration V stat computed at a constant speed and then to compute the vibrations around this configuration. The weak formulations of both sub-problems are nearly identical and relies on [Géradin and Rixen 1993, Curnier 1994, Batoz and Dhatt 1996, Ladevèze 1996[START_REF] Laîné | Dynamique des rotors[END_REF]]. Let us consider the transformation from the initial configuration V p0q into the current configuration, say V ptq (it is in particular valid for V stat ) as described in Fig. 3.2.

Since the in-service load is assumed known for the current configuration V ptq, the weak formulation is first written for this configuration then recovered to the undeformed configuration V p0q, which is commonly the one that is meshed (and machined and mounted in the engine). In all that follows, any quantity relative to the current deformed configuration is subscripted by x and any quantity relative to the undeformed configuration is subscripted by p. The initial configuration V p0q is supposed at rest, therefore the initial conditions are t 0 ùñ 6 9 9 8 9 9 7

x p, u p pp, 0q 0, (3.1) fu p ft pp, 0q 0.

For any configuration V , let ξpV q be the vector space of kinematically admissible displacement fields on this configuration. ξpV q is a subspace of the Sobolev space H 1 pV q (see Appendix C.2 for details on functional spaces).

The weak form of the mechanical problem relative to the current configuration

V ptq is § § § § §
Find u x px, tq ξpV ptqq such that dw x pxq ξpV ptqq, W a pu x , w x q W i pu x , w x q W e pu x , w x q,

where

W a pu x , w x q » V ptq ρ x pxq xa x px, tq , w x pxqy dV ptq (3.3)
is the virtual work of the acceleration quantities,

W i pu x , w x q ¡ » V ptq σ x pu x px, tqq : ε x pw x pxqq dV ptq (3.4)
is the virtual work of the interior forces and

W e pu x , w x q » V ptq
xf Vx px, tq , w x pxqy dV ptq » fV ptq

xf Sx px, tq , w x pxqy dSptq

(3.5)
is the virtual work of the exterior forces. σ x is the Cauchy strain tensor, ρ x is the volumic mass, a x is the acceleration, f Vx is the exterior volumic forces, f Sx is the exterior surfacic effort and x. , .y is the canonical inner product of ξpV q. All the latter quantities refer to the current configuration V ptq.

Remark: choice of the virtual displacement

The transformation that turns V p0q into V ptq between times 0 and t depends only on the real displacement u and not on w. Therefore, it is always possible to define a virtual displacement w such that w is

• either a function defined on the current configuration V ptq, that is a function of x. It is then denoted w x pxq and it belongs to H 1 pV ptqq, • or a function defined on the initial configuration V p0q, that is a function of p. It is denoted w p ppq and it belongs to H 1 pV p0qq.

There exists the following equivalence:

w p ppq w x pxq w x pxpp, tqq.

(3.6)

The time t has no importance for w since the transformation of w p into w x only requires a change of spatial coordinates, even if the coordinates depend on the time t at which the current configuration is chosen!

The aim of this section is to recover the terms in Eqns. (3.3) to (3.5) to the undeformed configuration V p0q so that the weak form can be written § § § § §

Find u p pp, tq ξpV p0qq such that dw p ppq ξpV p0qq, W a pu p , w p q W i pu p , w p q W e pu p , w p q .

(3.7)

The transformation of the different quantities that appear in Eqns. (3.3) to (3.5) are detailed in Secs. 3.1.2 to 3.1.4 respectively. Prior to that, a few definitions relative to transformations are recalled in Sec. 3.1.1 to help to clarify the notations.

General definitions

V ptq is the deformed structure and V p0q the same structure undeformed (see Fig. 3.2). The particle that is in p at time t 0 is in xpp, tq at time t, in other words x is the mapping that transforms V p0q into V ptq between times 0 and t. The initial condition is then xpp, 0q p.

The gradient of the transformation is denoted F p and is defined such that F p pu p pp, tq, tq fx fp pp, tq D p xpp, tq.

(

The notation

D v ¦ f fv ¦
where v is a field will be used throughout the rest of this chapter. If one defines the displacement u p pp, tq such that u p pp, tq xpp, tq ¡ p, F p pu p pp, tqq is also

F p pu p pp, tq, tq I D p u p pp, tq, (3.9)
where I is the idempotent mapping. The Jacobian of the transformation is denoted Jpu p pp, tq, tq and defined such that Jpu p pp, tq, tq det F p pu p pp, tq, tq.

The dilatation tensor (or right Cauchy-Green tensor) C p is defined by (3.12)

C p pu p pp
This tensor is referred to as the Cauchy's infinitesimal strain tensor or Green tensor.

Remark: small strain assumption

When the current configuration V ptq can be considered close to the initial one V p0q, the quantities defined previously become:

x p, F p pu p pp, tq, tq I, Jpu p pp, tq, tq 1, E p pu p pp, tq, tq ε p pu p pp, tq, tq.

( 3.13) This assumption is commonly used for vibration analysis.

Acceleration quantities: centrifugal/gyroscopic effects

The virtual work of the acceleration quantities W a pu x , w x q involves the acceleration a x . In the (fixed) reference frame, the velocity is given by [START_REF] Desceliers | Dynamique non linéaire en déplacements finis des structures tridimensionnelles viscoélastiques en rotation[END_REF]]

v x px, tq fx ft Ωptq x, (3.14)
or by v p pp, tq fu p ft pp, tq Ωptq pp u p pp, tqq.

(

This equation states states that the elastic motion described by the term fup ft pp, tq exists around a rigid body motion of rotation described by the term Ωptq pp u p pp, tqq, where the current position involves the elastic motion. The acceleration is derived from the previous expressions leading to (3.17)

a x px
Ωptq contains the information relative to the rotation of the structure. If ω x ptq, ω y ptq and ω z ptq are the components of the rotation vector ωptq relative to a basis pe x , e y , e z q, the matrix product

Ωptq¤ is equivalent to the vector product ωptq where Ωptq is the skew-symmetric matrix

Ωptq ! 0 ¡ω z ptq ω y ptq ω z ptq 0 ¡ω x ptq ¡ω y ptq ω x ptq 0 ( ) .
(3.18)

The expression of the acceleration in the reference frame exhibits different terms that are easy to interpret:

• the acceleration in the rotating frame f 2 up ft 2 pp, tq; • the centrifugal acceleration Ωptq 2 pp u p pp, tqq;

• the Coriolis acceleration dΩptq dt pp u p pp, tqq 2Ωptq fup ft pp, tq.

The virtual work of the acceleration quantities is finally recovered to the initial configuration (3.20)

V p0q W a pu p , w p q » V p0q
One can derive a set of bilinear and linear forms with respect to u p and w p . They are accompanied by the corresponding matrices and vectors derived from a Ritz-Galerkin procedure that relates the continuous field u p pp, tq to the discrete vector tqptqu via the matrix of shape functions rNppqs. These forms are: (i.) the bilinear form of mass, M pu p , w p q » V p0q ρ p ppq xu p pp, tq , w p ppqy dV p0q, A few remarks can be made after these definitions:

• for a given axis of rotation, clockwise and counterclockwise directions have opposite effects on the gyroscopic coupling, which is not the case for centrifugal softening. Centripetal and centrifugal accelerations have however opposite effects on the centrifugal acceleration.

• the inertial load is fully described by the sum K a pu p , w p q K c pu p , w p q F ei pw p q (resp. prK a ptqs rK c ptqsq tqu tf ei ptqu). It is decomposed into two distinct parts:

a follower part that changes as the state of deformation changes, with an explicit linear dependence on u p (resp. tqu). It is described by the bilinear forms of gyroscopic coupling, centrifugal acceleration and centrifugal softening.

a non-follower part that remains constant (modulo the changes in rotation speed) since it is computed on the reference configuration. It is contained in the global exterior force.

Remark: time-dependence of the bilinear and linear forms

All the bilinear and linear forms defined in Secs. 3.1.2 to 3.1.4 depend strongly on t in a fairly complex way. In particular, the forms that involve the rotation speed depend on t through the law imposed for its evolution in time (startup, shutdown, small variations around the nominal regime, and so on ...).

Exterior forces: follower forces

The virtual work of the exterior forces W e pu p , w p q is recovered to the initial configuration V p0q by where n p pp, tq is the outward-pointing normal of the elementary surface dSp0q.

W e pu p , w p q » V p0q xf Vp pp,
The exterior forces accompany the deformation of the structure, they are said to be follower forces. This is typically the case for pressure loads and inertial load.

Remark: inertial load as follower forces

The expression of the virtual work of the acceleration quantities involves volumic inertial terms that can be interpreted as inertial loads:

• the inertial centrifugal effort f ic ¡ρ p ppq Ωptq 2 pp u p pp, tqq ,

• the inertial gyroscopic effort f ig ¡2ρ p ppq Ωptq fu p ft pp, tq ¡ ρ p ppq dΩptq dt pp u p pp, tqq .

In the rotating frame, they are considered as exterior forces. One can see that these quantities follow the geometry, in other words their value is updated as the geometry changes. This is why they are often referred to as follower forces. Inertial load is often treated like this in finite element softwares [START_REF] Msc | MSC Software[END_REF]].

As a result, a linear form with respect to w p can be defined and attached to a vector representation in a subsequent Ritz-Galerkin procedure: (vi.) the linear form of exterior force, For the sake of brevity, the quantities relative to the global exterior force are now denoted F e pw p q and tf e ptqu and are given by F e pw p q F ext pw p q ¡ F ei pw p q, tf e ptqu tf ext ptqu ¡ tf ei ptqu .

F ext pw p q ¡ » V p0q
(3.34)

Interior forces: non-linear stiffness

The virtual work of the interior forces requires the recovery of σ x pu x px, tq, tq : ε x pw x pxqq to the initial configuration V p0q. Writing down the explicit calculation of this quantity leads to σ x pu x px, tq, tq : ε x pw x pxqq σ x pu x px, tq, tq : D x pw x pxqq t , since σ x is symmetric. Moreover, the derivative of w with respect to x is such that D x pw x pxqq D p w p ppq pF p pu p pp, tq, tqq ¡1 . This allows to write (3.39)

W i pu p , w p q ¡ » V p0q ¡ σ x pu x px,
A few more forms can be derived from Eqn. (3.39) to which matrix and vectors can be attached:

(vii.) the linear form of interior force which remains non-linear with respect to u p pp, tq, tq,

F i pu p , w p q » V p0q
rF p pu p pp, tq, tq S p pu p pp, tq, tqs : D p w p ppq dV p0q, (3.40) associated with the interior force vector,

tf i pqptqqu » V p0q
rD p Nppqs t : rF p pqptqqs rS p pqptqqs dV p0q. (3.41) The weak form of the interior forces can be developed further if one knows how to relate the stress tensor to the displacement, or at least to strain tensor. It involves the constitutive equation of a linear isotropic elastic material which is S p pu p pp, tq, tq H p ppq : E p pu p pp, tq, tq, (3.42) where H p ppq is the fourth-order tensor associated with Hooke's law of elasticity.

The substitution of S p by the constitutive equation in Eqn. (3.39) (3.45)

The time-dependence of G p is due to the fact that in visco-elasticity, stresses depends on the history of strains and inversely.

The weak problem described by Eqn. (3.7) is non-linear. It is equivalent to the research of the zeroes of a certain function R eq , called the residue in equilibrium, which turns the continuous problem to be solved into § § § § §

Find u p pp, tq ξpV q such that dw p ppq ξpV q, R eq pu p , w p q 0, (

where

R eq pu p , w p q M ¢ f 2 u p ft 2 , w p D g
¢ fu p ft , w p K a pu p , w p q K c pu p , w p q F i pu p , w p q ¡ F e pw p q .

(3.47) Equation (3.46) can then be replaced by its discrete form § § § § §

Find tqptqu C N such that d twu C N , tR eq pqptq, wqu t0u . (3.48) Since the residue R eq pqptq, w, tq is linear with respect to twu, Eqn. (3.48) can be turned into twu t tR eq pqptqqu t0u .

( 3.49) This equation must be satisfied whatever twu in C N . It leads to the following vector problem tR eq pqptq, tqu t0u .

(3.50)

The discrete form of Eqn. (3.47) returns tR eq pqptqqu tR eq pqptqqu rMs tX qptqu rD g ptqs t Wqptqu rK a ptqs tqptqu rK c ptqs tqptqu tf i pqptqqu ¡ tf e ptqu , (3.51) where t Wqptqu Both these non-linear minimization problems are usually solved with a linear iteration scheme (Newton-like algorithms). The idea is to reach the non-linear equilibrium by solving the linear problem tangent to the current configuration at each step. Let V 0 be a configuration characterized by u 0p (resp. tq 0 u) and δu p be an increment of the displacement around this state of equilibrium. The operator tangent to the residue R eq pu p , w p q at u 0p is defined by D up R eq u 0p ¨$ pδu p pp, tq, w p q K a pδu p , w p q K c pδu p , w p q D up F i pu 0p q $ pδu p pp, tq, w p q.

(3.52)

The matrix tangent to the residue vector tR eq pqptqqu at tq 0 u is defined the same way by rD q R eq pq 0 qs rK a ptqs rK c ptqs rK N L pq 0 qs .

(3.53)

The tangent operator is called (ix.) the bilinear form of tangent stiffness, (3.55)

D u0 p F i pu 0p q $ pδu p pp, tq, w p q » V p0q
The latter bilinear form (resp. matrix) states how the stiffness of the structure is affected by changes in its geometry. A remarkable property is that the limit when u 0p (resp. tq 0 u) tends towards 0 (resp. t0u) of the bilinear form of tangent stiffness (resp. the tangent stiffness matrix)

is the bilinear form of elastic stiffness

K e pu p , w p q » V p0q
D p w p ppq : H p ppq : D p u p pp, tq dV p0q. rD p Nppqs t : rH p ppqs : rD p Nppqs dV p0q).

(3.57)

Summary: quasi-static and vibration problems

The previous forms and matrices allow to compute both the quasi-static deformed configuration and the vibratory state around this configuration. Let the time t be momentarily suspended and suppose that one aims to determine the modal characteristics of the structure at an operating point where its angular speed Ωptq and angular acceleration W Ωptq are frozen at constant values now denoted Ω and W Ω. Matrices rD g s, rK a s, rK c s and vector rf e s are constant quantities at this operating point.

(i.) One aims to determine first how much the structure is statically deformed due to pΩ, W

Ωq. If any, other static loads such as thermal or gravity loads are also taken into account at this step. From Sec. 3.1.2, the centrifugal load is divided into follower and non-follower terms. The problem described by Eqn. (3.50) is first considered so that one can determine the static state described by the discrete displacement tq st u which nullifies the residue in quasi-statics tR eq pq st qu prK a s rK c sq tq st u tf i pq st qu ¡ tf e u . (3.58) tf e u contains only the static loads, in particular the non-following inertial load tf ei u whose expression is recalled here

tf ei u » V p0q ρ p ppq rNppqs t rΩ 2 W
Ωs tpu dV p0q. (3.59) This problem is solved with either a linear algorithm (in one step with rK e s) or a non-linear algorithm (iteratively with rK N L pq 0 qs updated at some steps, not necessarily at all steps)

until the configuration V stat is reached (see Fig. 3.1).

Remark: prestressed state of cyclically symmetric structures

When dealing with bladed disks or rotors, Chap. 1 highlights the convenience of the Fourier decomposition. With proper assumptions in the case of a rotor, this leads to a one-or two-sector per disk approach for either statics or dynamics.

The static-stress state is induced by the inertial load whose pattern is driven by δ 0, provided that the volumic mass is the same in all sectors of a given disk. Consequently, if the disk is tuned or if the rotor is made of tuned disks, the static displacement is sought as a static solution with δ 0.

Remark: time-dependent rotation speeds

The resolution techniques described in this section are also valid for a varying rotation speeds provided that a complementary discretization in time is made, for instance with the help from Newmark related schemes. At each time step t k , a quasi-static configuration induced by pΩpt k q, W Ωpt k qq is computed as the zero of a residue such that of Eqn. (3.58) with the help from a Newton-like algorithm. This provides the initial guess of the Newton algorithm at the following time step t k 1 . Such studies are beyond the scope of this work and are not considered here.

(ii.) Now that the quasi-static computations have returned the static state of deformation characterized by tq stat u, the vibratory response of the rotating structure has now to be determined around this static configuration. Vibrations are supposed to occur with small amplitudes so that the vibration problem normally described by Eqn. (3.50) is exactly the linear problem tangent to the non-linear problem at tq stat u. In this case, the discrete displacement tqptqu is the increment in displacement tδqptqu in Sec. 3.1.4 assumed small to comply with the vibration framework. The free or forced response satisfies the problem in Eqn. (3.50), in other words is sought as a zero of the operator tangent to the residue rD X q R eq pq st qs tX qptqu rD Wq R eq pq st qs t Wqptqu rD q R eq pq st qs tqptqu t0u .

(3.60)

The explicit matrix problem to be solved is then rMs tX qptqu rD g s t Wqptqu prK a s rK c s rK N L pq st qsq tqptqu tf ext ptqu , (3.61) where tf ext ptqu contains the exterior dynamic loads, since their static counterparts are balanced by the static solution. It is in particular the case of the non-follower part of the inertial load tf ei u. Note however that the follower part of this load is still present through rD g s, rK a s, rK c s and rK N L pq st qs. Finally, particular attention must be paid to the fact that tqptqu is said to depend on t, as Ωptq and W

Ωptq do, even if the latter are fixed for a modal analysis at a so-called operating point.

Parametrization in rotation

The axis of rotation is assumed to be e z . Note that it is always possible to define this axis such that The key idea of parametrization is to look how the different matrices involved in the static and dynamic problems are related to the rotation speed.

Linear matrices and vectors

The transformation of the matrices and vectors that are linear with respect to the rotation speed is straightforward. It leads to the collection of parametric matrices (3.64)

rD g s 2 Ω » V p0q

Polynomial stiffness matrix

In the following, one assumes that the rotation speed is constant (stabilized regime). Consequently W Ω 0, rK a s r0s and the second term in tf ei u is zero.

rK N L pq st qs and tf i pq st qu depend non-linearly on tq st u. In the general case, the dependence of tq st u on Ω is not obvious except if the static-stress state is close enough to the undeformed configuration so that one can neglect the follower forces, hence the effects of rK c s, and the non-linear part in tf i pq st qu. The last assumption is that the other static forces are negligible with respect to the non-follower part of the inertial load, so that only tf ei u remains. The quasi-static problem is then linear rK e s tq st u tf ei u . (3.65) It is solved in only one iteration and Eqn. (3.64) (3.68)

The coefficients in the polynomial given in Eqn. (3.67) are exact when tq st u satisfies Eqn. (3.66) and rK c K N L pq st qs is exactly a quadratic polynomial of the squared rotation speed. Otherwise, thanks to its simplicity, rK c K N L pq st qs is approximated by a quadratic matrix polynomial of the squared rotation speed Ω whatever its range.

For thin structures at low rotation speed or massive structures at moderate to high rotation speeds, it is commonly observed that the rotation speed changes the modal behaviour only little, hence it is common to perform computations with a static linear solver before the modal analysis (restart of a solution sequence 101 into a solution sequence 103 in MDNastran [START_REF] Msc | MSC Software[END_REF]). Consequently the polynomial shown in (3.67) is well suited to mode computations around the static-stress state. It is very convenient in terms of numerical cost since for parametric studies it requires only the assembly and storage of the matrix coefficients, which justifies the use of such an approximate form.

Remark: Other static loads

Actually, one should keep the other static loads contained in tf ext u, especially when one considers a constant thermal state (but not necessarily uniform). As a result, the expression of tq st u shown in Eqn. (3.66) is supplemented by a constant solution so that tq st u rK e s ¡1 ¡ tf ext u ¡

Ω 2 tf Ω ei u © tq Ω,0 st u Ω 2 tq Ω,1 st u, (3.69)
The matrix coefficients are then somewhat more complex and are not reproduced here for clarity.

Determination of the matrix coefficients

In other words, rK c K N L pq st qs is assembled exactly for given values of Ω then interpolated by a quadratic polynomial of the squared rotation speed between these points such as in Eqn. (3.67).

In the following, let rKpΩqs be rK c K N L pq st qs (the explicit dependence of rK c s and tq st u on Ω is omitted for clarity).

The determination of the matrix coefficients with Eqn. (3.68) is possible only when the assumption of proximity of the deformed state to the undeformed one holds. The matrix coefficients are to be determined as the solution of a matrix interpolation problem. The coefficients in the polynomial that interpolates each coefficient in matrix rK c K N L pq st qs are determined in a single operation. In general, if the highest degree in the matrix polynomial is N p , it has N p 1 matrix coefficients. The model has to be fully assembled (static prestress plus non-linear matrices) at N p 1 rotation speeds pΩ 0 , . . . , Ω Np q leading to the matrix equation dj v0, N p w, Np p0 Ω 2p j rK Ω,p s rKpΩ j qs (3.70) which is a collection of N 2 independent scalar equations (if rKpΩqs has a size of N ¢ N ). Using the Kronecker tensor product it can be written under the matrix form rΥ Is (3.73)

" " " ! K Ω,0 K Ω,1 . . . K Ω,Np ( 0 0 0 ) " " " ! KpΩ 0 q KpΩ 1 q . . . KpΩ Np q ( 0 0 0 ) , ( 3
A necessary condition for rΥs to be invertible is that its determinant is not zero. This latter K Ω,0 K Ω,1 . . . .76) In practice, the polynomial approximation of rK c K N L pq st qs must match the exact quadratic polynomial in Ω 2 at the limit when tq st u tends towards t0u, which implies that N p 2. The determinant of rΥs is consequently

condition is satisfied if dpi, jq N 2 Np , i $ j ùñ |Ω i | $ |Ω j |.
K Ω,Np ( 0 0 0 ) Υ ¡1 I $ " " " ! KpΩ 0 q KpΩ 1 q . . . KpΩ Np q ( 0 0 0 ) . ( 3
det prΥsq pΩ 2 2 ¡ Ω 2 1 q pΩ 2 2 ¡ Ω 2 0 q pΩ 2 1 ¡ Ω 2 0 q.
(3.77)

The matrix coefficients can be given under their analytical form

rK Ω,0 s

Ω 2 1 Ω 2 2 pΩ 2 2 ¡ Ω 2 0 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 0 qs ¡ Ω 2 0 Ω 2 2 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 1 qs Ω 2 0 Ω 2 1 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 2 ¡ Ω 2 0 q rKpΩ 2 qs , rK Ω,1 s ¡ Ω 2 1 Ω 2 2 pΩ 2 2 ¡ Ω 2 0 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 0 qs Ω 2 0 Ω 2 2 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 1 qs ¡ Ω 2 0 Ω 2 1 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 2 ¡ Ω 2 0 q rKpΩ 2 qs , rK Ω,2 s 1 pΩ 2 2 ¡ Ω 2 0 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 0 qs ¡ 1 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 1 ¡ Ω 2 0 q rKpΩ 1 qs 1 pΩ 2 2 ¡ Ω 2 1 q pΩ 2 2 ¡ Ω 2 0 q rKpΩ 2 qs . (3.78)
When the range of rotation speed is r0, Ω max s, with

Ω 0 0, Ω 1 1 2 Ω max , Ω 2 Ω max , the matrix coefficients are rK Ω,0 s rKp0qs , rK Ω,1 s 1 3Ω 2 max ¢ 16 K ¢ 1 2 Ω max & ¡ rK pΩ max qs ¡ 15 rKp0qs , rK Ω,2 s 4 3Ω 4 max ¢ rKpΩ max qs ¡ 4 K ¢ 1 2 Ω max & 3 rKp0qs . ( 3

.79)

As expected, the linear stiffness of the structure appears directly in the matrix polynomial as the term rK Ω,0 s.

Summary: prestressed vibration problem

From Secs. 3.1.5 and 3.2.3,the prestressed vibration problem of Eqn. (3.61) is turned into rMs tX qptqu Ω rD Ω g s t Wqptqu rK Ω,0 Ω 2 K Ω,1 Ω 4 K Ω,2 s tqptqu tf ext ptqu .

(3.80)

The matrices that appear in Eqn. (3.80) are constant matrices. The explicit dependence on the rotation speed allows to perform straightforward parametric computations.

Evolution of eigenfrequencies: Campbell diagrams

Parametric modal analyses are commonly considered for a rotating structure to determine how its eigenfrequencies and modeshapes are affected by the rotation speed. Many authors have studied the effects of the different matrices involved in such problems:

• effects of the gyroscopic coupling matrix rD g s are studied in [START_REF] Vidoli | Veering Phenomena in Systems With Gyroscopic Coupling[END_REF],

• effects of the centrifugal softening matrix rK c s are studied in [Chen and Peng 1995, Jacquet- Richardet et al. 1996, Huang and Kuang 2001, Moyroud et al. 2002, Tsai 2004, Chatelet et al. 2005],

• effects of the centrifugal acceleration matrix rK a s are studied in [Huang and Kuang 2001],

• effects of the non-linear tangent stiffness rK N L s are studied in ref. [Chen and Peng 1995, Tsai 2004, Petrov and Ewins 2005].

One is first interested in the stiffening/softening effects. rK c s tends to decrease the overall stiffness, this is particularly clear in Eqn. (3.64). Conversely, rK N L pq st qs tends to stiffen the overall stiffness, with the limit case of no additional stiffening when Ω, hence tq st u, is zero. The problem is that these two adverse effects act at the same time. In two early papers summarized in [Bazoune 2005], Southwell proposed that the eigenfrequencies of a rotating plain disk are bounded by two quadratic interpolations so that

ω 2 j pΩq ¥ ω 2 j p0q Ω 2 and ω 2 j pΩq ¤ ω 2 j p0q ω S,j Ω 2 , (3.81)
with ω S,j the Southwell coefficient associated with mode j. With such an approximation, the stiffening of the disk is the only possible behaviour because a given eigenfrequency may only increase with respect to the rotation speed. From this point of view, the matrix polynomial interpolation introduced in Sec. 3.2.3 produces more accurate results since it accounts for both softening and stiffening effects, because it reproduces the sum rK c K N L pq st qs which allows both types of evolution (increase or decrease). In fact, the eigenvalues are free to increase, decrease or switch between both in a given speed range.

The bladed disk of Fig. 2.1 serves as an illustration. It is assumed tuned so that complex conjugate pairs of mono-harmonic modes with Fourier coefficient δ 1 can be computed whatever the rotation speed thanks to Eqn. (3.80). The static state is computed as the response to an excitation with δ 0. Another phenomenon which has a strong influence is frequency veering. It refers to the phenomenon that sees two eigenvalues of a parameter-dependent system converge and then veer apart without crossing as this parameter evolves [START_REF] Afolabi | On Curve Veering and Flutter of Rotating Blades[END_REF]Mehmed 1994, Marugabandhu andGriffin 2003]. Veering is different from coalescence since the coupling of modes is conservative in the case of veering and non-conservative in the case of coalescence, the latter is typically involved in flutter [Balmès 1993, Afolabi and[START_REF] Afolabi | On Curve Veering and Flutter of Rotating Blades[END_REF]. Frequency veering is observed in Fig. 3.3b where double frequencies ω 6 and ω 7 (if they are numbered at Ω 0) associated with complex conjugate modeshapes tΦ 6 u, t s Φ 6 u ¨and tΦ 7 u, t s Φ 7 u come close around Ω 7.5 then veer away from each other. If the modeshapes are numbered at each speed point, modeshape tΦ 6 u (resp. t s Φ 6 u) before the veering becomes modeshape tΦ 7 u (resp. t s Φ 7 u) after the veering and inversely [START_REF] Balmès | High modal density, curve veering, localization: a different perspective on the structural response[END_REF]].

Extensive study of gyroscopic effects are beyond the scope of this work. rD g s is skewsymmetric and the resulting modes are complex. When a pair of modes exist at Ω 0, the gyroscopic coupling separates this pair into a forward-and backward-travelling wave as soon as Ω ¡ 0. Therefore, the Campbell diagram of such a structure displays two distinct curves as soon as Ω ¡ 0 [Jacquet- Richardet et al. 1996] instead of the two superimposed curves observed when the gyroscopic coupling is neglected, as it was done in the present section.

Remark: mistuning and rotation speed

When mistuning acts on the physical parameters (mass, stiffness, damping ...), it affects evenly rK a s, rK a s and rK e s, thus rK N L pq st qs.

The matrix polynomial is built from such matrices computed at three different speeds. For a tuned system without reduction, they are those of the initial sector. In the case of a mistuned system, these matrices are those of the whole disk model, reduced or not.

The mistuning pattern has to be applied to rKpΩqs, prior to the determination of the matrix coefficients. Anyway, as far as the author knows, detailed studies about how mistuning evolves with respect to the rotation speed are not available in the literature.

Multi-model approaches

Now that it has been demonstrated (resp. assumed) that the behaviour of a rotating structure relies exactly (resp. approximately) on a set of constant matrices, it is possible to use multi-model approaches.

Principle and key features

Dynamic problems such that described by Eqn. (3.80) form parametric families. They can be solved using a multi-model (MM) reduction technique [Balmès 1996b]. The key idea is to consider that their solutions are sought in an approximate subspace generated by target solutions computed at selected values of the parameter, here the rotation speed Ω. In other words, the solution tqptqu of Eqn. (3.80) at any speed Ω is sought as a linear combination of exact solutions tq j ptqu computed at N Ω target speeds Ω j , tqptqu N Ω j1 η j tq j ptqu . (3.82) tq j ptqu generally belongs to a subspace spanned by a basis rTpΩ j qs. Therefore, tqptqu is directly sought as a linear combination of the vectors in the union of the bases tqptqu rT Ω s tη Ω ptqu, (3.83) where tη Ω ptqu is a vector that contains the generalized DOF and

rT Ω s

T pΩ 1 q ¤ ¤ ¤ T pΩ n q $ (3.84)
spans a constant subspace. Columns in rTpΩ j qs, j v0, N Ω w, are target vectors that will be reproduced exactly by the reduced multi-model [START_REF] Balmès | High modal density, curve veering, localization: a different perspective on the structural response[END_REF]]. As far as the author knows, the first tentative to introduce a multi-model approach is found in [Marugabandhu and Griffin 2003] who used the subspace spanned by the modes of the structure at Ω 0 to generate the solution at any other rotation speed.

However, rT Ω s taken as is generates a subspace but may be made of linearly dependent vectors. Orthogonalization/elimination strategies similar to those introduced in Sec. 2.3.2 are employed to make the vectors linearly independent so that they really form a basis. Besides, the elimination process leads to a final set of generalized DOF not necessarily related to all target rotation speeds, even though all target vectors are later contained exactly within the subspace.

As a result, it seems difficult to exhibit an analytical expression of the generalized coordinates tη Ω u. Anyway, they are determined as the solution of a reduced matrix problem.

If one considers the structural problem in Eqn. (3.80), any finite element matrix rAs is projected onto the subspace generated by rT Ω s such that rA Ω s rT t Ω A T Ω s.

(3.85) This is in particular true for the stiffness matrix polynomial. As it is quadratic with respect to Ω 2 , N Ω 3 in Eqn. (3.84) is necessary to achieve sufficient accuracy, that is to say at least three target speeds have to be considered. This is very convenient since the determination of the matrix coefficients in the stiffness polynomial also requires to assemble matrices at three rotation speeds. It is possible, but not mandatory, to perform both the assembly of matrices and a subsequent modal analysis in a single operation at each target rotation speed.

When it comes to its applications, the multi-model approach can be employed (i.) on its own. When one is interested in the evolution of a set of mono-harmonic modes, the properties shown in Chap. • the different sets of vectors rTpΩ j qs, j v1, N Ω w, can be the outputs of distinct reduction processes performed at each target speed Ω j . The sets of vectors that span all these subspaces are then assembled together in rT Ω s before the vectors are made linearly independent. Section 3.3.4 presents an example that aims to build a multi-model Craig-Chang-Martinez reduction technique.

Exact Campbell diagrams of the HEBD model

The first illustration concerns the bladed disk of Fig. 2.1 considered as tuned. Thanks to Chap. 1, all modes are recovered from a family of parametric problems driven by the Fourier harmonic coefficient δ. The Campbell diagram relative to the modes of the disk is built by gathering the individual Campbell diagrams of the families of mono-harmonic modes. As a consequence:

• modes with different Fourier harmonic coefficients δ can cross but not veer,

• modes inside the same set of mono-harmonic modeshapes can veer but not cross. This is illustrated in Fig. 3.4 where the Campbell diagrams of modes with δ 1 and δ 2 are superimposed. These two families cross, but modes veer inside their own family, as displayed in Fig. 3.4b. As an example, the Campbell diagrams of the mono-harmonic modes with δ v0, 11w in the frequency band r0, 26s are computed within the speed range r0, 10s. The true operating range is r0, 1s, however it was stretched up to Ω 10 so that veering could be observed. The reference diagram for all harmonics is displayed in Fig. 3.5a to give an idea of the very high modal density in this frequency band (recall that modes come in pairs for δ ¡ 0). The first family of modes is selected and displayed in Fig. 3.5b, it corresponds to the first bending motion of the blade. A characteristic parabolic shape is observed at low speeds, which makes it compatible with classical observations [Marugabandhu and Griffin 2003]. For each individual harmonic coefficient δ v0, 11w: (i.) a parametric multi-model is built with a basis of modes computed at Ω t0, 5, 10u and Ω t0, 7.5, 10u. The number of modes retained for target speed equals the number of modes in the frequency band r0, 26s at Ω 0. The observed veering in speed [START_REF] Afolabi | On Curve Veering and Flutter of Rotating Blades[END_REF] does not affect the multi-model and modes are retrieved in the same order as the reference modes. As displayed in Fig. 3.7, the error stays low for mode pairs 6 and 7 that veer around Ω 7.5 and the modeshapes are perfectly correlated for every value of speed in general and around the veering points in particular. The multi-model and reference Campbell diagrams are very close to each other, as confirmed by the fact that the maximum relative error is either under 0.7% except for the last pair of modes with δ 5 in Fig. 3.8a or under 1.2% except for the 10 th pair of modes with δ 2 in Fig. 3.8b. This also illustrates the fact that the choice of the speed values to build the basis is the key point: Ω 7.5 as the intermediate speed increases accuracy. However, an acceptable default choice is

Ω 1 Ω min , Ω 2 1 2 pΩ min Ω max q and Ω 3 Ω max ,
it produces results that are accurate enough. Anyway, the error is zero for the target modes corresponding to the retained values of Ω in the basis.

Multi-model CICR

The computations in the previous section require the resolution of 3 ¢ 12 problems to build the multi-model basis with all harmonics as targets and cannot render mistuned configurations.

The results of Secs. 2.4 and 3.3.1 allow the use of a subset of target harmonics, that is to say a subset of mono-harmonic solutions with particular values of δ, to generate either the non-target tuned modes or the mistuned modes.

The multi-model CICR aims to consider a set of such target solutions computed at three rotation speeds as the input of a subsequent CICR so that it will generate

• in the tuned case, (i) the modes of the same target families at the non-target rotation speeds and (ii) the non-target tuned modes at any rotation speed,

• in the mistuned case, the modes at any rotation speed.

The mass matrix and the coefficients in the stiffness polynomial are projected onto the created subspace and parametric reduced models are solved to compute the Campbell diagrams and the forced responses in both tuned and mistuned cases.

As an illustration, δ v0, 3w are defined as target harmonic coefficients and Ω t0, 5, 10u

as target speeds, with the same constraint as in 3.3.2 for the number of retained modes. The Campbell diagram of the tuned system is very close to Fig. 3.5a and will not be reproduced here for brevity. The maximum relative error on frequency is plotted for each speed and each value of δ. The target modes are indeed found exactly by the multi-model as displayed in Fig. 3.9a and the error stays below 0.02% in the frequency band of the target modes. Figure 3.9b also shows that the maximum of the relative error committed on frequency is still small, below 1.4%, in a three times wider frequency band. Two distinct sources of error appear on this graphs. The first error is committed in the generation of the non-target tuned modes, it can be seen in the planes Ω 0, Ω 0.5 and Ω 1. Note that the error in plane Ω 0 in Fig. 3.9b is exactly the enveloppe of the plot in Fig. 2.14. The second error, characterized by typical bell curves in Fig. 3.9a, is committed in the generation of all the modes at the non-target speeds. in the model to ensure that the Campbell diagram for the modes with δ 1 can be quite well superimposed to the response in Fig. 3.10a (a loss factor η 10 ¡3 is chosen arbitrarily). As expected, the tuned response only involves the family of modes with the same Fourier harmonic coefficient. Figure 3.10b clearly shows that mistuning affects some modes by splitting their double frequencies. Some regions are particularly affected by frequency scattering and avoiding resonances in these regions may be a particularly difficult task.

Multi-model Craig-Chang-Martinez reduction

The second example is related to a booster blade (denoted BB in the following) whose root is assumed clamped. It is displayed in Fig. 3.11. The study aims to build a reduced multi-model condensed at three nodes located at the tip of the blade for subsequent contact analyses (which are beyond the scope of the present study). Two steps are required: (i.) at each target speed a reduction basis is built so that the model is condensed at the chosen master DOF by the means of a Craig-Chang-Martinez reduction technique;

(ii.) the three sets of vectors are orthogonalized so that the multi-model reduction basis is obtained.

The first step consists in the condensation of the model at the nine retained DOF. This type of condensation is generally achieved by a Craig-Bampton reduction technique. In this framework, the kinematic subspace is approximated by

• a set of constraint modes, which are static responses to successive unitary deflections of each interface DOF,

• fixed interface modes.

This corresponds to a reduction basis that is rTs

I 0 ¡K ¡1 cc K ci Φ f ix & . ( 3

.86)

Fixed interface modes lead to a loss of information at the interface in dynamics. In particular, a non negligible error can be committed on normal modes of a model condensed to the DOF set of the interface when the latter is free. To address this issue, the latter subspace is enriched by substituting fixed interface modes by free interface modes. This variant is the Craig-Chang-Martinez technique. The approximation subspace is rTs

I ¡K ¡1 cc K ci Φ & . ( 3

.87)

This set of vectors is linearly dependent due to the fact that it contains too much information, therefore it is not full rank. It is orthogonalized so that it becomes rTs

I 0 ¡K ¡1 cc K ci Φ c K ¡1 cc K ci Φ i & . ( 3

.88)

In the present multi-model approach, one has to build one Craig-Chang-Martinez per operating speed due to the presence of the stiffness matrix in the reduction basis. To lighten notations, one defines

rK cci s K ¡1 cc K ci $ , rΦ Ω s Φp0q Φ 1 2 Ω max ¨Φ pΩ max q $ .
(3.89)

For further orthogonalization purposes, one gathers, sorts and divides the vector sets

rT i s I ¡K cci p0q & , rT c s I I ¡K cci 1 2 Ω max ¨¡K cci pΩ max q Φ Ω & . ( 3

.90)

A difficulty arises because the generalized DOF of the constraint modes in each set are related to the physical DOF of the interface as stated by the identity matrices. In other words, the constraint modes at each speed refer to the same set of generalized DOF. It requires to choose

¡K cci 1 2 Ω max ¨ K cci p0q ¡K cci pΩ max q K cci p0q Φ Ωc K cci p0q Φ Ωi & . (3.91)
Nevertheless, vectors within rT c s are not themselves linearly independent. They are processed through an iterative maximum sequence (see [START_REF] Balmès | Orthogonal Maximum Sequence Sensor Placements Algorithms for modal tests, expansion and visibility[END_REF]] and Sec. 2.3.2 for details) to become so. The final multi-model basis is returned by gathering rT i s and rT c s in rT Ω s

T i T c $ . ( 3

.92)

This reduction technique is applied to determine the modal characteristics of the booster blade within the speed range r0, 1s. Figure 3.12a displays the Campbell diagram used as reference.

The true operating range of this blade is r0, 0.25s and it was designed to present a fairly flat Campbell diagram in this region. Anyway, the range of operating speed is virtually extended up

to Ω 1 so that the influence of the speed on the eigenfrequencies is high enough to give rise to veering events. The correlation between the modeshapes at Ω 0 and Ω 0.5 (resp. Ω 1)

indicates that the eigenmodes are quite affected in the chosen range.

Table 3.1 provides a comparison in execution time when at each iteration in speed:

• a problem built with the multi-model Craig-Chang-Martinez reduction technique with Ω t0, 0.5, 1u as the target speeds is solved (MMCCM), • a problem built with the multi-model Craig-Chang-Martinez reduction technique with Ω t0u as the target speed is solved (CCM0),

• the full matrix problem is assembled then solved (Full).

In the third case, the matrix polynomial is used to build the stiffness matrix at the current operating point. The two multi-model Craig-Chang-Martinez reduction techniques seems to produce results in a significant smaller amount of time. Their accuracy is directly related to the contents of the initial set of vectors. As depicted in Fig. 3.13a, the results are excellent in the case of a basis computed with Ω t0, 0.5, 1u since the error stays below 2.5% for the first 30 modes, the worst case being achieved for the two modes that are the closest to the upper bound of the frequency band (modes 29 and 30). Figure 3.13a displays a characteristic bell curve, stating that, as expected, the error is forced at zero for the target modes at Ω 0, Ω 0.5 and Ω 1, and the error is greater between these values. In the second case, when the basis is computed with Ω t0u only, the multi-model diverges rapidly from the reference model as can bee seen in Fig. 3.13b where 50% of error is reached for Ω 0.14. Moreover, numerical problems are noticed as soon as Ω reaches 0.2.

One is still concerned about the stability of the proposed approach to frequency veering. One considers the veering event that occurs between modes 22 and 23 near Ω 0.81 (see Fig. 3.14a).

The error shows that, as expected, the multi-model is not affected by veering as it returns the modes in the same order as the reference.

Conclusions

This chapter first details the formulation employed to deal with structures rotating at a constant angular speed Ω. The finite element model derived from this weak formulation exhibits additional gyroscopic and stiffening/softening terms that affect the free and forced response of the structure in both their frequencies (eigenfrequencies or resonances) and shapes. The numerical cost to assemble the full problem at each desired rotation speed can be prohibitive for structures clearly affected by inertial load. A first approximation is proposed where the stiffness matrix is turned into a quadratic polynomial of the squared rotation speed. This polynomial is exact in the limit case when the statically deformed structure is close to the undeformed one. This operation simplifies the assembly of the stiffness matrix but keeps the size of the subsequent matrix problem unchanged. Nevertheless, it allows straightforward parametric studies of the vibration characteristics of rotating structures.

The stiffness polynomial can be used in a multi-model approach that helps to reduce the latter size. It assumes that the solution in displacement can be found in a kinematic subspace generated by solutions computed exactly for a set of three target speeds. As a result, the size of the problem is reduced down to the size of the target set of vectors once they have been made linearly independent. The accuracy is controlled by the size of this set. The multi-model approach is very useful to draw precise Campbell diagrams since it becomes rapidly profitable. Another advantage is its stability to veering events. Moreover, it is obviously compatible with other reduction techniques. When used as an initial step, it only affects the input set of vectors and only the coefficients in the matrix polynomial have to be projected onto the reduced subspace. Practical uses of this methodology are illustrated on various examples.

Methodological aspects

The flowchart in Fig. 3.15 summarizes the multi-model reduction technique. A variant exists where the polynomial is assembled prior to the computations of the target modes. The polynomial is evaluated at each target speed and a subsequent modal analysis is performed to retrieve the mode set.

Further developments

The present study paves the way for extensive parametric studies with the help from multimodel approaches. What is presented in this chapter could be extended by

• including gyroscopic coupling so that stability studies can be conducted. Gyroscopic coupling depends explicitly on the rotation speed (see Sec. 3.1.4) and multi-model approaches are thus straightforward.

• keeping an eye on promising concepts introduced in the field of non-linear dynamics such as Non-Linear Normal Modes (NLNM), which are defined as (not necessarily synchronous) periodic solutions of a non-linear system [START_REF] Peeters | [END_REF]]. Instead of first searching for the non-linear quasi-static equilibrium of the structure under the inertial load and then computing its eigenmodes, one could be interested in computing directly its NLNM. One of the particularities of these modes is the fact that their frequency and shape depend on the amount of energy injected to the structure. In the present case, a fraction of such energy comes from the inertial load, proportionally to the squared rotation speed. This would be a straightforward determination of the vibration characteristics of the rotating structure. However, for the moment, NLNM cannot be used as is for modal synthesis which currently limits their practical use.

• extending the multi-model approach to other parameters such as the temperature in the structure or the pressure applied to the blades. The temperature is a critical issue for rotors since there exists a thermal gradient globally from the turbine to the front stage of compressor and locally from the blade to the inner rim of the disk.

Conclusions and Perspectives

The general outcome of this study is the introduction of methods allowing the prediction of multi-stage full rotor dynamics on the basis of refined 3D models currently used for detailed stress analysis. These methods are demonstrated to be practically applicable to current industrial models.

The first contribution is a review of the Fourier decompositions of spatially periodic geometries and its extension to the independent transformation of stages in the case of full rotors. The matrix form of the Discrete Fourier Transform is shown to lead to the classical decomposition into independent sub-problems for a so-called tuned stage (a model whose mechanical properties are identical from one bladed sector to another). When the latter assumption fails, full finite element matrices of the mistuned system are involved in the matrix problem because of inter-harmonic coupling. Particular coupling patterns are presented and identified according to the distribution of the mechanical properties in the assembled matrices. Current rotor designs lead to modes that involve inter-stage coupling and full rotor predictions are thus required. The assembly of geometrically periodic structures that belong to different groups of symmetry is however generally not periodic itself. By considering that disks are meshed independently and linked by an intermediate ring that has no particular condition of symmetry, the DFT can be applied at the disk level and the inter-harmonic coupling is shown to be entirely linked to the energy in the inter-stage coupling ring. It is also shown that coupling through a continuity condition generates similar problems. The concept of mono-harmonic multi-stage modes [Laxalde et al. 2007a] therefore appears as an approximation where the off-diagonal inter-harmonic coupling terms are neglected. Detailed examples are used to demonstrate that in practice the true system modes are often very close to being mono-harmonic (the approximation is generally valid), some cases of multi-harmonic behavior are however described. Current work on multi-stage analysis of rotors with visco-elastic damping devices using constrained layers with low number segmentation will be an opportunity to further analyze the aliasing issue found in mono-harmonic modes. Indeed, mono-harmonic solutions are shown to be subject to aliasing, with high diameter modes that do not converge towards the fixed disk blade frequencies. The analysis of the spatial content of full rotor modes (made possible by the methods introduced in the second contribution) demonstrates that this aliasing is actually a distortion with the true modes continuing to exhibit the classical single stage accumulation of frequencies. This problem illustrates that the exact synthesis of a full rotor response from the combination of spatially mono-harmonic responses may not be possible in the case of multi-stage assemblies.

The second contribution is the introduction of reduction techniques that use target modes to create reduced models of the full disk. The target modes are classical cyclic symmetry solutions in the case of single stage assemblies and mono-harmonic solutions in the case of full rotors. Blade modes, where the sector edges are fixed, are also introduced to guarantee the proper prediction of accumulation frequencies for high diameter modes. The fundamental interest of the approach is that based upon the computation of a few cyclic solutions (typically less than 5 diameters), one creates a fairly compact model of the full rotor. It is a reduced model and can therefore be used to predict modes, free and forced responses, etc. . The reduction bases enable full and partial restitution and this has been optimized for very large models. The model predicts the retained diameter modes exactly and gives a very good approximation of all the other. The resulting model is also a rotor model with superelements representing sectors individually. Starting from a nominally tuned model, it is thus possible to replace individual sectors by models with other properties. The usefulness of this ability is demonstrated by a disk mistuning analysis, where one considers that the variable parameter is the blade frequency of individual sectors. Probabilistic analyses would be a straightforward extension and multi-stage mistuned rotor analyses would certainly be a useful application. It should also be noted that the approach is perfectly compatible with a direct replacement of the model of a given sector and could thus easily be used to study "rogue" blades such as found after a foreign object damage for instance. The idea of building a reduction basis from target modes is first developed by ensuring compatible displacements on the sector left and right interfaces. This approach is shown to be useful for single disks and to allow mistuning analyzes. Its extension to multi-stage assemblies would however lead to a large number of DOF to ensure continuity. A second approach is therefore introduced that enforces coupling by reducing disjoint components (no common degree of freedom) linked by interfaces that have a physical extent (typically a single layer of elements between the components). In the implementation one distinguishes inter-sector interfaces, which are repeated and thus represented as superelements, and inter-disk interfaces which are only used once and conserved as standard elements. Collinearity is a fundamental difficulty when building reduction bases from arbitrary sets of vectors. Various orthogonalization methods are thus discussed. In particular, one extends the iterative maximum sequence method proposed in [START_REF] Balmès | Orthogonal Maximum Sequence Sensor Placements Algorithms for modal tests, expansion and visibility[END_REF]] to allow a sorting of the DOF sets according to their coupling nature and generations of bases with particular patterns of zero contributions on target DOF sets. The various procedures considered have highlighted that a trade-off exists between basis size and final sparsity in the model. The need to preserve certain features exactly, such as blade modes, has also been seen. Refinement of the orthogonalization strategies is certainly needed to guarantee the quality of predictions while obtaining the smallest possible reduced rotor model. The scalability of the procedure is illustrated on an industrial full 3D rotor model with 65 million DOF for which no direct computations are possible. A significant effort was placed on illustrating procedures that can be used to analyze the full rotor response. Direct analyses of generalized responses include indicators such as the contribution of each sector of each disk to the overall response, and transforms to view the spatial harmonic content. The second indicator is in particular used to test whether the true responses are really mono-harmonic or not. Recovery of the response to viewing meshes or partial meshes for energy or stress analyses are also possible. The analysis of stresses was limited to simple verifications of dynamic stress distribution in certain modes. The procedures used to compute multi-stage static responses to centrifugal and thermal loads, which are necessary for fatigue analyses (drawing of Haigh diagrams), have been implemented but are not reported here. Partial recovery on viewing meshes can be used to analyze modes and responses but also generated forced responses to point or distributed pressure loads. This has been used to validate a test/analysis correlation strategy for modal analysis results available at Snecma. The finalization of this correlation is a short-term perspective.

The third contribution of this work is the use of the multi-model approach [Balmès 1996b] to generate reduction bases that accurately predict modes for a range of rotation speeds at a low numerical cost. It is shown that the geometric stiffness matrix and inertial softening effects can be represented by a quadratic polynomial of the squared rotation speed. This result is exact when the deformed geometry is close enough to the initial one and an approximation otherwise. Using target modes that are exact solutions computed at three different speeds, one can build a reduced model that gives very accurate results at all intermediate speeds. This process is shown to be totally adapted to mono-disk and multi-stage methods introduced earlier: the input set of vectors normally computed at a single speed is simply complemented with solutions computed at two other speeds, with no impact on the implementation of the orthogonalization procedures. Illustrations demonstrate the accuracy including the ability to represent veering events. A straightforward extension would be the inclusion of gyroscopic effects leading to differentiation of modes into rotating waves. The issue would then be to validate the need to further enrich the bases to account for the gyroscopic loads which are proportional to the rotation speed. Another possibility would be to consider thermal loads as a parametric variation, the interpolation of matrices by the use of discrete thermal states does not however appear as straightforward as the use of a polynomial of the rotation speed, so that the usefulness of the idea remains unclear. when the Laplace variable s is s i χ.

2. pχ, p ypχqq is called the spectrum of function y, it is continuous. 3. In signal processing, x is the time variable, then χ is called the (natural) frequency. When dealing with geometrically repetitive structures (by translation or rotation), x is either a linear or angular position and χ is then either a linear or angular frequency. More generally, if y belongs to C n pRq, n N ¦ and v is such that v y pnq , p vpχq i n χ p ypχq.

Proof 1.

1. The proof is trivial.

2. The proof of the first property requires an integration by parts, but here no term yp0q arises since the boundaries of the integral are ¡V and V. Then, the formula for the derivative of order n is obtained by a recurrence scheme initiated by this property.

3. The discrete Fourier transform is also virtually periodic with respect to χ with a period of N s ∆χ or equivalently the sampling rate f s . Indeed, The spectrum is generally represented for χ varying between 0 and N s ∆χ, but it can be equivalently plotted for χ varying between ¡N s ∆χ{2 and N s ∆χ.

Sampling complex-valued signals whose frequencies are of the form δ ∆χ m N s ∆χ δ ∆χ m f s , m Z, at sample rate f s leads to identical sets y s , n v0, N s ¡ 1w, hence to the same spectrum p y δ , k v0, N s ¡ 1w. Such signals are said aliases of each other. When they are real-valued, signals whose frequencies are of the form ¨k ∆χ m N s ∆χ, m Z, are aliases with respect to f s . In other words, a given line at k 0 ∆χ with 0 ¤ k 0 N s {2 may correspond to:

(i.) a signal whose frequency is χ 0 k 0 ∆χ, (ii.) a signal whose frequency is χ 0 pk 0 m N s q ∆χ, (iii.) a signal whose frequency is χ 0 pm N s ¡ k 0 q ∆χ.

As a result, to avoid the presence of aliases, the sampled signal must be such that its maximum frequency χ max is below half the sampling rate f s :

χ max f s 2 , (B.8)
or on the contrary, the sampling rate must be at least twice the maximum frequency of the sampled signal if the latter is guessed. This theorem is known as the Shannon-Nyquist theorem and f s {2 is commonly referred to as the Nyquist frequency.

This phenomenon is illustrated in

B.1.3 Real-valued sampled functions

In this section, y is assumed real-valued, i.e. such that y : R ÝÑ R

x Þ ÝÑ ypxq and the series y s is real. In the following for convenience, the direct and inverse Discrete Fourier transforms are rewritten separating real and imaginary quantities as Repp y δ q cos ps δ αq ¡ Impp y δ q sin ps δ αq .

( (B.27)

It can be recovered under a matrix form by tyu t tyu tp yu t rE t Es tp yu tp yu t tp yu.

(B.28)

B.1.5 Generalized orthogonality conditions

The following orthogonality conditions are not needed for the Fourier transform. However The condition of definiteness does not hold since for any u H 1 pV q such that dx V , fu fx pxq 0 does not imply that dx V , upxq 0, because any constant function on V satisfies this property. Therefore, for any v H 1 pV q, xu , vy K,V 0 whereas u $ 0.

Theorem C.2.5. K-semi-norm

||.|| K,V : H 1 pV q ÝÑ F u Þ ÝÑ ||u|| K,V § § § § § § § § fu fx § § § § § § § § M,V £ » V ¢ fu fx pxq 2 dV 1 2
.

is a semi-norm.

Proof 5. The condition of separation does not hold either for the same reason. Therefore, it is possible that ||u|| K,V 0 whereas u $ 0.

Remark 11.

1.

¡ H 1 pV q, x. , .y M,V © is a Hilbert space.

2.

¡ H 1 pV qztu, dx V , upxq a, a F 3 zp0, 0, 0qu, x. , .y K,V © is a Hilbert space, since the form

x. , .y K,V is an inner product on this space only.

Figure 1 :

 1 Figure 1: A few gas turbines for various uses: (a) high bypass ratio turbofan, (b) low bypass ratio turbofan, (c) turboprop, (d) helicopter turbine and (e) ground turbine

Figure 2 :

 2 Figure 2: (a) A high bypass turbofan (GP7200 -© Engine Alliance) and (b) variation of the parameters across the engine

Figure 5 :

 5 Figure 5: Haigh diagram with the approximate limit curves and safe areas (positive stresses only) -σ S : static stress, σ D : dynamic stress, σ Y : yield stress, σ U : ultimate stress and N c : number of cycles

Contents 1 . 1 Fourier

 11 Decomposition on a single periodic geometry . . . . . . . . . 13 1.1.1 Fourier transform of fields . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.2 Fourier transform of discretized fields . . . . . . . . . . . . . . . . . . . 15 1.1.3 Fourier transform of matrices . . . . . . . . . . . . . . . . . . . . . . . . 16 1.1.4 Continuity conditions between sectors . . . . . . . . . . . . . . . . . . . 17 1.1.5 Response of a tuned disk to external excitations . . . . . . . . . . . . . 20 1.1.6 Constraint elimination . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.1.7 Norms of the discretized fields . . . . . . . . . . . . . . . . . . . . . . . 25 1.1.8 Comments on the methodological aspects . . . . . . . . . . . . . . . . . 26 1.2 Mistuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.2.1 Characteristics of mistuning . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.2.2 Phenomenological aspects . . . . . . . . . . . . . . . . . . . . . . . . . . 30 1.3 Multi-stage dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.3.1 Assembly of two disks through an intermediate ring . . . . . . . . . . . 35 1.3.2 Hypothesis of decoupled harmonics . . . . . . . . . . . . . . . . . . . . . 39 1.3.3 Response of a rotor made of tuned disks to external excitations . . . . . 41 1.3.4 Assembly of two disks with a constraint equation . . . . . . . . . . . . . 43

Figure 1 . 1 :

 11 Figure 1.1: Sample bladed disk

Figure 1 . 2 :

 12 Figure 1.2: Lumped model of a single tuned disk

Figure 1 . 3 :

 13 Figure 1.3: Projected matrices of a single tuned lumped disk with 12 sectors

Figure 1

 1 Figure 1.4: Sub-domains of sector S s

  .21) (ii.) more generally, if the continuity condition described by Eqn.(1.16) is assembled on the full disk, this results in

  1.1.

Figure 1 . 5 :

 15 Figure 1.5: Eigenfrequencies vs. Fourier harmonic coefficients for a sample disk with 23 blades: ( ) disk-dominated motion and ( ) blade-dominated motion with (-) the corresponding accumulation frequencies and modeshapes

Figure 1

 1 Figure 1.6: Response of a sample disk with 23 blades to an excitation with δ 2: (a) full response with (--) observed resonances and (b) zoom to peaks 1M with (¤ ¤ ¤ ) the frequencies of all the modes in that band

Figure 1 . 7 :

 17 Figure 1.7: Projected matrices rMs and rKs of a single tuned lumped disk with 12 sectors after elimination of the inter-sector continuity conditions

Figure 1 .

 1 Figure 1.10: Tuned and mistuned responses of a sample disk with 23 blades to an excitation with δ 2 and DFT of the response at some peaks when the distribution has a standard deviation of (a) 0.5%, (b) 5%, (c) 10% and (d) 20%

  Similar results are obtained with a cosine mistuning pattern cosps ζ αq thanks to Eqn. (B.29).

Figure 1 .

 1 Figure 1.11: Projected stiffness matrix r p K ∆ p Ks of a lumped disk with 12 sectors when (a) tuned and with (b) sine harmonic mistuning with ζ 1, (c) alternate mistuning (cosine harmonic mistuning with ζ 6) and (d) random mistuning.

Figure 1 .

 1 Figure 1.12: History of multi-stage computations at Snecma: (a) Fourier axisymmetric, (b) 3D slice, (c) 3D cyclic symmetry, (d) 360°shell/volume model and (e) multi-stage cyclic symmetry

Figure 1 .

 1 Figure 1.13: Sample model of a multi-stage assembly

Figure 1 .Figure 1 .

 11 Figure 1.14: Lumped model of an assembly of tuned disks

Figure 1 .Figure 1 .

 11 Figure 1.16: Projected matrices of an assembly of disks with 5 and 8 sectors respectively: (a) matrices of disk D 1 , (b) matrices of disk D 2 and (c) matrices of ring R

Figure 1 .

 1 Figure 1.17: Harmonic correspondence pattern for a rotor with (a) 4 and 8 sectors and (b) 5 and 8 sectors: (-) disk D 1 and (-) disk D 2

Fig. 1 .

 1 Fig.1.18b, since lcmp5, 8q 0 and due to the fact that the harmonic correspondence pattern described in Fig.1.17 is much longer, many more blocks are filled. Aliasing issues are discussed further in Sec.1.3.3.

Figure 1 .

 1 Figure 1.18: Harmonic correspondence patterns in r p K R s for a rotor made of tuned disks (a) with 4 and 8 sectors and (b) with 5 and 8 sectors (in the latter case, only the off-diagonal blocks are considered for legibility)

Figure 1 .

 1 Figure 1.19: Frequencies vs. Fourier harmonic coefficients for a rotor with 23 and 45 blades: (a) full graph with accumulations, ( ) multi-stage zone and ( ) aliasing zone and (b) zoom on f r0, 8.5s with the eigenfrequencies of the individual disks for various boundary conditions

  Let us consider two bladed disks D 1 and D 2 that are bonded together. The meshes of the aft rim I 1 a of D 1 and the front rim I 2 f of D 2 are generally different. Disks D 1 and D 2 have N 1 and N 2 DOF respectively. These two DOF sets are disjoint. Any matrix rAs relative to this rotor is given by Eqn. (1.56) in which rA R s is removed. Equation (1.58) without rA R s returns the projected matrix, which is block-diagonal and, from Sec. 1.3.3, does not contribute to intra-disk harmonic coupling. In the case of bonded disks, a constraint equation is required to reproduce: (i.) inter-sector continuity in each disk defined in Sec. 1.1.5. For each disk D d , d t1, 2u, it requires the definition of rc d r s, rc d l s and rc d s that appear in Eqns. (1.22) and (1.23), now superscripted by d.

Figure 1 .

 1 Figure 1.20: Lumped model of an assembly of bonded tuned disks

Figure 1 .

 1 Figure 1.21: Kernel basis of the constrained subspace of lumped rotors with (a) 4 and 8 sectors and (b) 5 and 8 sectors respectively

Figure 2 . 1 :

 21 Figure 2.1: Finite element model of an integrally bladed disk [HEBD]

Figure 2 . 2 :

 22 Figure 2.2: Academic integrally bladed rotor [ACAD]

Figure 2 . 3 :

 23 Figure 2.3: Industrial model of a whole rotor: (a) front and (b) rear [FSHP]

Figure 2

 2 Figure 2.4: Details of the rings (a) between D 1 and D 2 , (b) between D 2 and D 3 and (c) between D 3 and D 4 [FSHP]

Figure 2 .

 2 Figure 2.5: Sector models and rings used to compute the mono-harmonic eigensolutions (a) front and (b) rear [FSHP]

Figure 2 . 6 :

 26 Figure 2.6: Sample sector and disk models: (a) standard finite element model and (b) superelement model with (G l) interface DOF

Figure 2

 2 Figure 2.7: Iterative Maximum Sequence: (a) step 0, (b) step 1 and (c) step 2

Figure 2 . 8 :

 28 Figure 2.8: Typical selections for post-processing of (a) motion and (b) stresses.

Figure 2 .

 2 Figure 2.9: Basis [CICR]

Figure 2 .

 2 Figure 2.10: Reduced matrices rM s S s and rK s S s [HEBD]

Figure 2 .

 2 Figure 2.11: Reduced matrices rM S s and rK S s [HEBD]

Figure 2 .

 2 Figure 2.12: Pseudo-mesh of the generalized degrees of freedom [HEBD] According to Eqn. (2.20), when mistuning is small enough to keep rT s S,c s unchanged, it acts only on the diagonal terms rΛ s S,c s of the reduced stiffness matrix, in other words it corresponds to deviations of the natural frequencies of each sector S s with fixed interfaces, such that ds v0, N s ¡ 1w,

Figure 2 .

 2 Figure 2.13: Vectors in the reduction basis for the HEBD model: vectors in (a) rT s S,r s, (b) rT s S,c s and (c) rT s S,l s

Figure 2 .

 2 Figure 2.14: CICR technique: (a) Eigenfrequencies vs. Fourier harmonic coefficients and (b) relative error in frequency for the HEBD model

Figure 2 .

 2 Figure 2.15: Reduced representations of the modeshapes (left) and their recovery to the full HEBD (right) for mono-harmonic modes with: (a) δ 0, (b) δ 2 and (c) δ 11

Figure 2 .

 2 Figure 2.16: Fourier harmonic contents of the modes displayed in Fig. 2.15 with: (a) δ 0, (b) δ 2 and (c) δ 11

Figure 2 .

 2 Figure 2.17: Natural frequencies of the sectors with fixed interfaces at most ¨5% of the mean (tuned) value [HEBD]

  2.18 and 2.19 where two generalized modeshapes.

Figure 2 .

 2 18 depicts a mode mostly localized to blade #17, which is indicated equivalently by both representations. The reduced representation

Figure 2 .

 2 Figure 2.18: (a) Reduced representation of a localized modeshape of the HEBD, (b) its recovery to the full disk and (c) its harmonic content

Figure 2 .

 2 Figure 2.19: (a) Reduced representation of a degenerate tuned modeshape of the HEBD, (b) its recovery to the full disk and (c) its harmonic content

Figure 2 .

 2 Figure 2.20: Mistuned response of the HEBD to an excitation with δ 3 in the vicinity of pairs of tuned peaks number (a) 1, (b) 2, (c) 3, (d) 5, (e) 6, 7, 8 and (f) 9

Figure 2 .

 2 Figure 2.21: Sector substructuring in the mono-disk DJCR: (l) sector superelement, ( l) intersector superelement

Figure 2 .Figure 2 .

 22 Figure2.22 displays the transpose of a sample mono-disk reduction basis without vector ordering (actual basis) and with vector ordering (not used but shown here for clarity). The two triangular areas of non-zero data in Fig.2.22b correspond to typical results of the IMS with 1 as the diagonal terms.With such a topology, generalized DOF tq s 1 M,c u are not coupled with DOF of sector s since their motion on the right interface is zero. The coupling matrix rA s M,C s is thus significantly

Figure 2 .

 2 Figure 2.23: Sector substructuring in the multi-disk DJCR: (l) sector superelement, ( l) intersector superelement, ( l) ring surface and (G l) nodes common to the ring and the right interface

  s, non-zero shapes that have their maximum in the DOF set common to the right edge (I d,s r ) and the ring (R). These are in very limited number and thus kept non-zero everywhere. They are selected by a sorting of the initial vectors according to their contribution to the r R DOF set followed by a simple application of the IMS algorithm. • r r T d,s M,rR s, shapes that are only zero on the right edge/ring intersection r R. They are obtained by a simple application of the IMS algorithm to the remaining vectors. • other shapes are possibly non zero on both the interior S d,s c and ring DOF sets. To limit coupling, one builds r r T d,s M,R,R s by the singular value decomposition of the trace on R of the remaining vectors and r r T d,s M,r s by application of the IMS to the trace of the same vectors on c l. This step thus involves some replication, which was deemed acceptable.

Figure 2 .

 2 Figure 2.24: Sample multi-stage DJCR basis of a single sector with vectors (a) unsorted and (b) sorted [ACAD]

Figure 2 .

 2 Figure2.24 displays the transpose of a sample multi-disk reduction basis for disk D 1 without vector ordering (actual basis) and with vector ordering (not used but shown here for legibility). The triangular areas of non-zero data in Fig.2.24b correspond to typical results of the IMS.

Figure 2 .

 2 Figure 2.25: Topology of the reduced matrices according to the orthogonalization and DOF sorting strategy: (a) without and (b) with isolation of the DOF common to the ring and right interface -ring without mass (ACAD model with 3 and 4 sectors)

Figure 2 .

 2 Figure 2.26: Reduced stiffness matrix rK M s (a) global shape: ( ) blocks of disk D 1 , ( ) blocks of disk D 2 and ( ) blocks of ring R, and (b) actual matrix (ACAD model with 12 and 15 sectors)

Figure 2 .

 2 Figure 2.27: (a) Target frequency bands, (b) relative error on frequency and (c) modeshape correlation between reference and generalized modes (ACAD model with 12 and 15 sectors)

  2.28. Similarly to Fig. 1.19b in Sec. 1.3, one compares the rotor mono-harmonic eigensolutions to the eigenmodes of the individual disks in free-fix (D 1 ), fix-fix (D 2 ) and fix-free (D 3 ) conditions. Frequencies are normalized with respect to the frequency of the first flexible 0 diameter solution of the rotor. The same conclusions presented in Sec. 1.3 can be drawn from that graph and are briefly recalled: (a) accumulations are well reproduced except that (b) aliased solutions are observed for disks D 1 and D 2 for δ ¡ 18 and δ ¡ 28 respectively and (c) multi-stage coupling affects only the left-hand triangular area which concentrates the disk-dominated solutions (see Fig. 2.28b).

Figure 2 .

 2 Figure 2.28: (a) Full and (b) zoomed normalized frequencies of the mono-harmonic eigensolutions of the 3SHP model ( ) and the modes of disk D 1 ( ), disk D 2 (¢) and disk D 3 ( ¢) in free-fix, fix-fix and fix-free condition -(-), (-) and (-) are accumulations

Figure 2 .Figure 2 .Figure 2 .

 222 Figure 2.29: Coupled mono-harmonic eigensolution with δ 1 at f 15.53: (a) amplitude of the displacement and (b) density of strain energy -Blade shapes are 2S1, 2T and 2F [3SHP]

Figure 2 .

 2 Figure 2.32: Mono-harmonic eigensolution localized to disk D 3 with δ 0 at f 5.35: (a) amplitude of the displacement and (b) density of strain energy -Blade shape is 1F [3SHP]

Figure 2 .

 2 Figure 2.33: Reduced stiffness matrix rK M s of the 3SHP model (a) global shape: ( ) blocks of disk D 1 , ( ) blocks of disk D 2 , ( ) blocks of disk D 3 , ( ) blocks of ring R 12 , ( ) blocks of ring R 23 and (b) actual matrix

Figure 2 .

 2 Figure 2.34: Generalized modes 7 to 300 (f r1, 7.52s): (a) Spatial harmonic content and (b) participation factor in strain energy for disks ( ) D 1 , ( ) D 2 and ( ) D 3 [3SHP]

Figure 2 .

 2 Figure 2.35: Generalized modes 301 to 600 (f r7.52, 13.74s): (a) Spatial harmonic content and (b) participation factor in strain energy for disks ( ) D 1 , ( ) D 2 and ( ) D 3 [3SHP]

Figure 2 .

 2 Figure 2.36: Generalized modes 601 to 900 (f r13.76, 19.57s): (a) Spatial harmonic content and (b) participation factor in strain energy for disks ( ) D 1 , ( ) D 2 and ( ) D 3 [3SHP]

Figure 2 .

 2 Figure 2.37: Restored mode at f 10.23 localized to disk D 2 : (a) modeshape and (b) spatial harmonic content for disks ( ) D 1 , ( ) D 2 and ( ) D 3 [3SHP]

Figure 2 .

 2 Figure 2.38: Restored mode at f 15.89 with strong coupling: (a) modeshape and (b) spatial harmonic content for disks ( ) D 1 , ( ) D 2 and ( ) D 3 [3SHP]

Figure 2 .

 2 Figure 2.39: Typical selections for post-processing of (a) motion and (b) stresses [3SHP]

Figure 2 .

 2 Figure 2.40: Partial recovery of modeshapes at (a) f 10.23 and (b) f 15.89 [3SHP]

Figure 2 .

 2 Figure 2.41: Partial post-processing of the generalized mode at f 15.89 on disk D 1 : (a) normalized density of strain energy and (b) normalized distribution of Von Mises stress [3SHP]

Figure 2 .

 2 Figure 2.42: Normalized frequencies of the mono-harmonic eigensolutions of the FSHP model ( ) -(-), (-) and (-) are accumulation frequencies

Figure 2 .

 2 Figure 2.43: Global bending with δ 1 at f 1.57 [FSHP]

Figure 2 .

 2 Figure 2.44: Global torsion with δ 0 at f 1.69 -Blade shapes are 1F, 1F and 1F [FSHP]

Figure 2 .

 2 Figure 2.45: Coupled mono-harmonic eigensolution with δ 2 at f 13.45 -Blade shapes are 2T, 2F and 1T [FSHP]

Figure 2 .

 2 Figure 2.46: Mono-harmonic eigensolution localized to the seal with δ 1 at f 16.33 [FSHP]

Figure 2 .

 2 Figure 2.47: Reduced stiffness matrix rK M s of the FSHP model (a) global shape: ( ) blocks of disk D 1 , ( ) blocks of disk D 2 , ( ) blocks of disk D 3 , ( ) blocks of disk D 4 , ( ) blocks of ring R 12 , ( ) blocks of ring R 23 and ( ) blocks of ring R 34 and (b) actual matrix 1 867 modes lie in the normalized frequency band r0, 29s. They are computed in 7 854 s

Figure 2 .

 2 Figure 2.48: Evolution of the eigenfrequencies of the reduced model in the band r0, 24.17s with the accumulation frequencies

Figure 2 .

 2 Figure 2.49: Fraction of strain energy by spatial harmonic for the generalized modes in the frequency band r2.90, 6.44s relative to (from top to bottom) D 1 , D 2 , D 3 and D 4 [FSHP]

Figure 2 .

 2 Figure 2.50: Fraction of strain energy by spatial harmonic for the generalized modes in the frequency band r2.90, 6.44s relative to (from top to bottom for each harmonic) D 1 , D 2 , D 3 and D 4 [FSHP]

Figure 2 .Figure 2 .Figure 2 .Figure 2 .Figure 2 .Figure 2 .

 222222 Figure 2.52: Global bending at f 1.59: (a) partially recovered modeshape and (b) spatial harmonic content [FSHP]

Figure 2 .

 2 Figure 2.58: Graph f pδq with both initial and reconstructed solutions [FSHP]

  2.60. Recall that with the definition of the DFT given in Appendix B.1.3, harmonics between 1 and N d s {2 ¡ 1 also account for those between N d s {2 1 and N d s ¡ 1 (see Eqn. (B.12)), hence they seem to have a "double" amplitude.

Figure 2 .

 2 Figure 2.60: Harmonic content of the first single DOF load applied to the FSHP model

Figure 2 .

 2 Figure 2.61: Accumulation of modes in the vicinity of the (--) 1F frequency of disk D 3 at f 5.25 [FSHP]

Figure 2 .

 2 Figure 2.62: Effective stiffness ratios of the generalized modes for the first single DOF load: (¤ ¤ ¤ ) assumed controllable modes and (--) 1F accumulation frequency of disk D 3 [FSHP]

Figure 2 .

 2 Figure 2.63: (a) Maximum response on the blade DOF set of the pseudo-mesh and (b) MMIF for (-) disk D 1 , and (-) disk D 2 and (-) disk D 3 : (¤ ¤ ¤ ) assumed resonating modes and (--) 1F accumulation frequency of disk D 3 [FSHP]

Figure 2 .Figure 2 .

 22 Figure 2.64: Locations of the sensors used to perform a Fourier analysis to the transfer functions [FSHP]

Figure 2 .

 2 Figure 2.66: Measured displacement at sensors: (a) considered resonances and amplitudes at resonance (b) 3 and (c) 7 for (-) disk D 1 , and (-) disk D 2 and (-) disk D 3 [FSHP]

Figure 2 .Figure 2 .

 22 Figure 2.67: (a) Shape and (b) fraction of strain energy by harmonic at peak number 3 at f 3.99 [FSHP]

Figure 2 .

 2 Figure 2.69: Flowchart of a typical reduced mono-/multi-stage analysis within SDT

Figure 2 .

 2 Figure 2.70: Mistuned disk models made of superelements: (a) small mistuning and (b) geometric mistuning -(G l) interface DOF, (G l) interior DOF of the nominal sectors and (G l) interior DOF of the rogue sector

Contents 3 . 1

 31 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 3.1.1 General definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 3.1.2 Acceleration quantities: centrifugal/gyroscopic effects . . . . . . . . . . 115 3.1.3 Exterior forces: follower forces . . . . . . . . . . . . . . . . . . . . . . . 117 3.1.4 Interior forces: non-linear stiffness . . . . . . . . . . . . . . . . . . . . . 118 3.1.5 Summary: quasi-static and vibration problems . . . . . . . . . . . . . . 121 3.2 Parametrization in rotation . . . . . . . . . . . . . . . . . . . . . . . . 123 3.2.1 Linear matrices and vectors . . . . . . . . . . . . . . . . . . . . . . . . . 123 3.2.2 Polynomial stiffness matrix . . . . . . . . . . . . . . . . . . . . . . . . . 123 3.2.3 Determination of the matrix coefficients . . . . . . . . . . . . . . . . . . 125 3.2.4 Summary: prestressed vibration problem . . . . . . . . . . . . . . . . . 127 3.2.5 Evolution of eigenfrequencies: Campbell diagrams . . . . . . . . . . . . 127 3.3 Multi-model approaches . . . . . . . . . . . . . . . . . . . . . . . . . . 128 3.3.1 Principle and key features . . . . . . . . . . . . . . . . . . . . . . . . . . 129 3.3.2 Exact Campbell diagrams of the HEBD model . . . . . . . . . . . . . . 130 3.3.3 Multi-model CICR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 3.3.4 Multi-model Craig-Chang-Martinez reduction . . . . . . . . . . . . . . . 134 3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 3.4.1 Methodological aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 3.4.2 Further developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

Figure 3 . 1 :

 31 Figure 3.1: The two steps used to compute the response of a structure rotating at a constant speed

Figure 3 . 2 :

 32 Figure 3.2: Transformation of the initial configuration into the current one

F

  mass matrix, which is symmetric, definite and positive rMs » V p0q ρ p ppq rNppqs t rNppqs dV p0q, (3.22) (ii.) the bilinear form of gyroscopic coupling, D g pu p , w p q 2 » V p0q ρ p ppq xΩptq u p pp, tq , w p ppqy dV p0q, (3.23) associated with the gyroscopic coupling matrix, which is skew-symmetric rD g ptqs » V p0q 2 ρ p ppq rNppqs t rΩptqs rNppqs dV p0q, (3.24) (iii.) the bilinear form of centrifugal acceleration, K a pu p , w p q » ) the bilinear form of centrifugal softening, K c pu p , w p q » V p0q ρ p ppq e Ωptq 2 u p pp, tq , w p ppq i dV p0q, (3.27) associated with the centrifugal softening matrix, which is symmetric, definite and positive rK c ptqs » V p0q ρ p ppq rNppqs t rΩptqs 2 rNppqs dV p0q, (3.28) (v.) the linear form of non-follower exterior inertial load,

ρ

  p ppq xf Vp pp, tq , w p ppqy dV p0q ¡ » fV p0q xf Sp pp, tq , w p ppqy dSp0q, (3.32) associated with the exterior force vector, tf ext ptqu ¡ » V p0q ρ p ppq rNppqs t rNppqs tf V p ptqu dV p0q ¡ » fV p0q rNppqs t rNppqs tf Sp ptqu dSp0q.(3.33)

  the linear elastic stiffness matrix which is symmetric, semi-definite and positive rK e s » V p0q

Figure 3 . 3

 33 displays a typical Campbell diagram of this family of modes. Their true behaviour is quite far from the quadratic curves they would describe if they were interpolated as in Eqn. (3.81).

Figure 3 . 3 :

 33 Figure 3.3: Campbell diagram of mono-harmonic modes with δ 1 (a) speed range r0, 10s and (b) zoom on the veering region [HEBD]

Figure 3 . 4 :

 34 Figure 3.4: Campbell diagram of mono-harmonic modes with δ t1, 2u (a) speed range r0, 10s and (b) zoom on a veering region [HEBD]

Figure 3

 3 Figure 3.5: Campbell diagram of mono-harmonic modes for δ v0, w11 (a) for all modes (with engine orders 1 to 11 in the background) and (b) for the 1F family of modes [HEBD]

Figure 3 .

 3 6 displays an excellent correlation in frequencies and modeshapes between the resulting Campbell diagram and the reference one for modes with δ 1 when the multi-model basis is built at Ω t0, 5, 10u. Note that the Modal Assurance Criterion invoked here is that of Eqn. (1.47) since the modeshapes come in complex conjugate pairs.

Figure 3

 3 Figure 3.6: (a) Campbell diagram for δ 1 and (b) relative error on frequency for a multi-model basis built at Ω t0, 5, 10u -The color bar refers to the modeshape correlation [HEBD]

Figure 3

 3 Figure 3.7: (a) Campbell diagram for δ 1 and (b) relative error on frequency restricted to mode pairs 6 and 7 for a multi-model built at Ω t0, 5, 10u -The color bar refers to the modeshape correlation [HEBD]

Figure 3 . 8 :

 38 Figure 3.8: Maximum error on frequency for a multi-model basis at (a) Ω t0, 5, 10u and (b) Ω t0, 7.5, 10u [HEBD]

Figure 3 Figure 3 .

 33 Figure 3.9: Maximum error on frequency for a multi-model CICR of the HEBD model in the frequency band (a) r0, 5s and (b) r0, 15s

Figure 3 .

 3 Figure 3.11: Geometry of the booster blade

Figure 3 .

 3 Figure 3.12: Modes of the BB model (a) Campbell diagram and (b) MAC comparison of modes at Ω 0.5 and Ω 1 with modes at Ω 0.

Figure 3 .Figure 3 .

 33 Figure 3.13: Relative error on frequency and modeshape correlation for a multi-model Craig-Chang-Martinez reduction techniques built at (a) Ω t0, 0.5, 10u and (b) Ω t0u -The color bars refer to the modeshape correlation [BB]

Figure 3 .

 3 Figure 3.15: Flowchart of a typical reduced multi-model approach within SDT (with possible calls to MDNastran)

Fieldsε

  Linearized strain tensor σ Cauchy stress tensor τ Boussinesq stress tensor S Piola-Kirchhoff stress tensor (second) E Green-Lagrange strain tensor H Hooke tensor Operators θ Rotation of angle α D ¦ .

Remark 1 . 1 .

 11 Sets of numbersva, bw all integers between a and b ra, bs all real numbers between a and b ta, b, ¤ ¤ ¤ , zu set of integers/real numbers Fourier Transform Let y be a function defined byy : R ÝÑ C x Þ ÝÑ ypxq,such that y is integrable on R (in the sense of the distributions). The Fourier transform p y of y is defined The Fourier transform is a particular case of the bilateral Laplace transform,

1 . 2 .

 12 Linearity: given two functions py, vq integrable on R and pa, bq C 2 such that u a y b v, Fourier Transform of a derivative: given a function y in C 1 pRq and if v is such that v y I , there is the relation p vpχq i χ p ypχq.

  Fig. B.1 where three differents sines whose fundamental frequencies are χ 0 3, χ 0 7 and χ 0 13 respectively are sampled with a rate of f s 10 (hence ∆χ 1). With this sampling rate, the sines are aliases of each other since 3 3 ¢ ∆χ 0 ¢ f s , 7 ¡3 ¢ ∆χ 1 ¢ f s and 13 3 ¢ ∆χ 1 ¢ f s . Figure B.1b clearly shows that their spectra are undistinguishable.

Figure B. 1 :

 1 Figure B.1: Sampled sines with (-) χ 0 3, (-) χ 0 7 and (-) χ 0 13 at f s 10. (a) ypxq and (G) y s (b) |p y δ |.

Repp y δ q 1

 1

  Unitary inverse Fourier transform ds v0, N s ¡ 1w, q cos ps δ αq ¡ Impp y δ q sin ps δ αqq Notation are kept unchanged for the sake of legibility, but the latter definitions are those that are used throughout the main chapters of this manuscript. The unitary Fourier matrix rEs derives from Eqn. (B.25) as ppN s ¡ 1q δ αq ¡ c 2 sin ppN s ¡ 1q δ αq ¤ ¤ ¤ p¡1q Ns¡1 Parseval theorem For a real-valued sampled function, the use of unitary transforms lead to Parseval's theorem under q| 2 |Impp y δ q| 2 ¨.

ζζζ

  they are involved in the transformation of mistuned matrices. Let pδ, δ I , ζq N 3 . There are the relations: Ns¡1 ş0 4 cosps δ αq cosps ζ αq sinps δ I αq 0 Ns¡1 ş0 4 sinps δ αq cosps ζ αq cosps δ I αq 0 Ns¡1 ş0 4 cosps δ αq cosps ζ αq cosps δI αq $ 0 rN s s or N s ¡ δ δ I ζ rN s s or δ N s ¡ δ I ζ rN s s or 2 N s ¡ δ ¡ δ I ζ rN s s , 0 else.Ns¡1 ş0 4 sinps δ αq cosps ζ αq sinps δ I αq $ 0 δ $ 0 and δ $ N s {2 and δ I $ 0 and δ I $ N s {2 and 5 N s ¡ δ δ I ζ rN s s or δ N s ¡ δ I ζ rN s s , sinps ζ αq cosps δ I αq 0 Ns¡1 ş0 sinps δ αq sinps ζ αq sinps δ I αq 0 Ns¡1 ş0 cosps δ αq sinps ζ αq sinps δ I αq $ 0 δ I $ 0 and δ I $ N s {2 rN s s or N s ¡ δ δ I ζ rN s s or δ N s ¡ δ I ζ rN s s or 2 N s ¡ δ ¡ δ I ζ rN s s , sinps ζ αq cosps δ I αq $ 0 rN s s or N s ¡ δ δ I ζ rN s s or δ N s ¡ δ I ζ rN s s or 2 N s ¡ δ ¡ δ I ζ rN s s ,

r

  Zs tqu t f u,

	(1.36)
	with
	rT Ker s span pkerprcsqq , tqu rT Ker s tqu , r Zs T t $ , Ker Z T Ker t f u rT Ker s t tfu .

  if rAs rIs, since the identity matrix is definite, x. , .y I is an inner product and ||.|| I is a norm, called the Euclidean norm (or canonical norm). Notations x. , .y 2 and ||.|| 2 are

	2
	K 2 tqu t rKs tqu . 1	(1.42)
	(iii.) preferably employed.

  The inner product of two vectors tqu and tpu of R N is xq , py A tqu t rAs tqu tp qu rE I N 0 s t I Ns A 0 $ rE I N 0 s tp qu .

		(1.44)
	Thanks to Eqn. (B.22), the latter equation becomes
	xq , py A xp q 0 , p p 0 y A 0	Ns 2 ¡1 δ1

  , N s ¡ 1w, r∆A s s 2 pRe pr∆A ζ sq cosps ζ αq ¡ Im pr∆A ζ sq sinps ζ αqq , Im pr∆A ζ sq sinppN s ¡ 1q ζ αq $ rsin ζ ∆A ζ s rE I N 0 s rE t sin ζ Es r∆A ζ s. Sin ζ Es couple pairs of harmonics pδ, δ I q that satisfy one of

						(1.50)
	where r∆A ζ s is the generating matrix. With a single sine pattern, r∆As is
	r∆As ¡2	" " " " " " !	0 ¤ ¤ ¤ . . . . . . 0 ¤ ¤ ¤ Im pr∆A ζ sq sinps ζ αq ¤ ¤ ¤ 0 ¤ ¤ ¤ . . . . . . . . . . . . 0 ¤ ¤ ¤ 0 ¤ ¤ ¤ ( 0 . . . 0 . . . 0 0 0 0 0 0 )	,	(1.51)
	rsin ζ ∆A ζ s.
	This allows to compute r∆ p As which is given by
			r∆ p As	E t I N 0		(1.52)
	Equation (B.30) states that rE t the following conditions

Table 2 . 1 :

 21 Characteristics of the finite element model of the HEBD

	Stage	Nb. of	Nb. of	Nb. of	Nb. of	Type of
		sectors	Nodes	Elements	DOF	Elements
	1	1{23	4 785	809	14 355	Penta15 Hexa20
	1	23{23	97 819	18 607	293 457	Penta15 Hexa20

Table 2 . 2 :

 22 Characteristics of the finite element model of the academic rotor[ACAD]

	Stage	Nb. of	Nb. of	Nb. of	Nb. of	Type of
		sectors	Nodes	Elements	DOF	Elements
				Disks
	1	1{12	184		84	552	Hexa8
	1	12{12	1 872	1 008	5 616	Hexa8
	2	1{15	184		84	552	Hexa8
	2	15{15	2 340	1 260	7 020	Hexa8
	1 Ø 2	-	405	Ring	810	1 215	Tetra4

Table 2 . 3 :

 23 2.3. (ii.) the full rotor made of second order elements, referred to as the FSHP model. In this case, the finite element model is made of 12 456 113 elements (15 node pentahedrons, 20 Characteristics of the finite element model of the 3SHP

	Stage	Nb. of	Nb. of	Nb. of	Nb. of	Type of
		sectors	Nodes	Elements	DOF	Elements
				Disks
	1	1{36	16 343	65 170	49 029	Tetra4
	1	36{36	560 880	2 346 120	1 682 640	Tetra4
	2	1{56	9 566	36 815	28 698	Tetra4
	2	56{56	502 992	2 061 640	1 508 976	Tetra4
	3	1{70	23 104	97 860	69 312	Tetra4
	3	70{70	1 565 130	6 850 200	4 695 390	Tetra4
				Rings
	1 Ø 2	-	5 664 2 592	14 141 4 836	16 992 7 776	Tetra4 Tetra4
	2 Ø 3	-	3 094 2 562	7 515 4 782	9 282 7 686	Tetra4 Tetra4
	Table 2.4: Characteristics of the finite element model of the FSHP
	Stage	Nb. of	Nb. of	Nb. of	Nb. of	Type of
		sectors	Nodes	Elements	DOF	Elements
				Disks
	1	1{36	108 375	65 170	325 125	Tetra10
	1	36{36	3 803 148	2 346 120	11 409 444	Tetra10
	2	1{56	62 421	36 815	187 263	Tetra10
	2	56{56	3 382 176	2 061 640	10 146 528	Tetra10
	3	1{70	151 032	97 860	446 622	Tetra10
	3	70{70	10 387 090	6 850 200	31 161 270	Tetra10
	4 Ñ T	1{50	119 894	19 874	359 688	Penta15 Hexa20
	4 Ñ T 50{50	4 373 950	993 700	13 121 850	Penta15 Hexa20
				Rings
	1 Ø 2	-	21 264 8 640	56 853 20 735	63 792 25 920	Tetra4 Tetra4
	2 Ø 3	-	11 004 8 540	28 563 20 496	33 012 25 620	Tetra4 Tetra4
	3 Ø 4	-	15 400 4 600	65 210 12 596	46 200 13 800	Tetra4 Tetra4

node hexahedrons and 10 node tetrahedrons for the disks, and 4 node tetrahedrons for the ring) that connect 21 946 364 nodes, which leads to 65 839 092 DOF. The individual characteristics of the sectors and the stages are reported in Tab. 2.4.

As described in Sec. 2.2.2, one uses inter-disk rings built from a Delaunay tessellation between the node sets of the regarding rims of two adjacent disks. To ensure that the intermediate rings have no interior DOF, they are made of the first-order elements that are returned by the

  Ring R between disks D d and D d 1 connects tq d,s rR u and tq d,s R u, ds v0,N d s ¡ 1w, to tq d 1 ,s 1 rR u and tq d 1 ,s 1 R u, ds I v0, N d 1 s ¡ 1w.The restriction of the whole kinematic subspace to the DOF set of

	.	(2.34)
	the ring is

Table 2 .

 2 6: Memory occupation of the reduced matrices of the 3SHP model with 134 vectors in

	rT s init s
	Common DOF isolated	No	Yes
	Size	8 881	12 313
	Non-factorized	38.2Mb	50.5Mb
	Cholesky factor	628.9Mb 802.6Mb

Table 2

 2

	.7: CPU time needed to compute the mono-harmonic eigensolutions of the 3SHP model
	(s: single, np: near-pairs)
	δ	time (s)	modes	DOF
	0	258	62 s	140 565
	1	904	62 np	281 130
	2	847	58 np	281 130

  2.6.2. Figures 2.34, 2.35 and 2.36 in which the color varies from white (relative amplitude of harmonic equals 0) to black (relative amplitude of harmonic equals 1) then show that, as expected, modes come first by near-pairs, and second by families. The fraction of strain energy in each disk is defined by

Table 2 .

 2

	8: CPU time needed to compute the mono-harmonic eigensolutions of the FSHP model
	(s: single, np: near-pairs)
	δ	time (s)	modes	DOF
	0	29 590	79 s	1 318 692
	1	147 582	76 np	2 637 384
	2	125 621	68 np	2 637 384

Table 2 .

 2 9: Vector distribution in the initial set (s: single, np: near-pairs)[FSHP]

	Mono-harmonic eigensolutions	Fixed interface modes
	δ 0	79 s	D 1	9 s
	δ 1	76 np	D 2	5 s
	δ 2	68 np	D 3	3 s
	δ 3	58 np
	δ 4	47 np

Table 2 .

 2

		10: Characteristics of the reduction bases [FSHP]
	SE in Disk D 1 D 2 D 3 D 4	Basis size (intra ¢ reduced) 308 727 ¢ 184 174 156 ¢ 167 430 401 ¢ 166 228 378 ¢ 389	Number of interior DOF 71 58 45 110	Memory used 433.4 Mb 221.9 Mb 545.1 Mb 677.8 Mb

Assembling the reduced matrices requires 6 071 s, leading to matrices of size 47 668 ¢ 47 668, each one needs 3.8 Gb to be stored. Their typical topology is displayed in Fig.

2

.47. The different blocks are detailed in that figure.

Table 2 .

 2 11 summarizes the common and distinct features of the two proposed reduction techniques.

Table 2

 2

	.11: Common and distinct features of the two proposed reduction techniques
	CICR	DJCR
	Reduction type
	CMS with sector superelements	CMS with sector superelements
	Initial set of vectors
	Mono-harmonic modes	Mono-harmonic eigensolutions
	Modes of sector with inter-sector inter-	Modes of sector with inter-sector and
	faces fixed	inter-disk interfaces fixed
	Basis construction
	Orthonormalization of vectors via Gram-Schmidt	Sorting of vectors according to their con-tribution to each DOF set and orthogonal-ization DOF set by DOF set (SVD or IMS)
	Quite heavy vector handling since inter-sector continuity must be ensured	Light vector handling
		Generalized DOF set is related to the
	Generalized DOF are arbitrary	physical DOF set for straightforward assem-
		bly
	Multi-stage capabilities
		Yes, taken into account from the initial
	No	set of vectors (multi-stage mono-harmonic
		eigensolutions)
	Mistuning capabilities
	Yes, either small or large (the latter is	Yes, designed for either small or large
	not yet implemented)	but not yet implemented
	Multi-model capabilities
	Yes	Yes
	Enhanced post-processing
	Yes, including a direct visualisation of the generalized modeshapes	Yes

  , tq, tq F p pu p pp, tq, tq t F p pu p pp, tq, tq, C p pu p pp, tq, tq I D p u p pp, tq D p u p pp, tq t D p u p pp, tq t D p u p pp, tq. pu p pp, tq, tq t F p pu p pp, tq, tq ¡ I u p pp, tq D p u p pp, tq t D p u p pp, tq t D p u p pp, tq

						(3.10)
	From the latter tensor, one introduces the Green-Lagrange tensor E p such that
	E p pu p pp, tq, tq 1 2 2 E p pu p pp, tq, tq 1	©	,	©	.	(3.11)
	E

¡

F p ¡ D p p is symmetric and can be decomposed into a linear and a non-linear part. The linear part is denoted ε p and is ε p pu p pp, tq, tq

1

2 pD p u p pp, tq D p u p pp, tq t q.

  , tq fv x pp, tq f 2 u p ft 2 pp, tq dΩptq dt pp u p pp, tqq 2Ωptq fu p ft pp, tq Ωptq 2 pp u p pp, tqq.

	ft px, tq dΩptq dt	x 2Ωptq fx ft px, tq Ωptq 2 x,	(3.16)
	or to
	a p

  Jpu p pp, tqq ρ x pxq. Equation (3.19) can be expanded so that W a pu p , w p q

	» » V p0q V p0q » V p0q » V p0q	ρ p ppq 2ρ p ppq f f 2 u p ft 2 pp, tq , w p ppq f Ωptq fu p ft pp, tq , w p ppq p dV p0q p dV p0q ρ p ppq f¢ Ωptq 2 dΩptq u p pp, tq , w p ppq dt ρ p ppq f¢ Ωptq 2 dΩptq dt p , w p ppq	p	dV p0q

ρ p ppq xa p pp, tq , w p ppqy dV p0q, (3.19) since ρ p ppq p dV p0q.

  tq , w p ppqy dV p0q Vp pp, tq Jpu p pp, tq, tq f Vx px, tq, f Sp pp, tq Jpu p pp, tq, tq § § § § § §pF p pu p pp, tq, tqq ¡t n p pp, tq § § § § § § f Sx px, tq,

	»	fV p0q	xf Sp pp, tq , w p ppqy dSp0q, (3.31)
	if one defines
	f

  tq, tqpF p pu p pp, tq, tqq ¡t © : D p w p ppq Jpu p pp, tq, tq dV p0q ,(3.35) since dV ptq Jpu p pp, tq, tq dV p0q. From the expression above, one can derive the first Piola-Kirchhoff (or Boussinesq) stress tensor τ p defined by τ p pu p pp, tqq Jpu p pp, tq, tq σ x pu x px, tq, tq pF p pu p pp, tq, tqq ¡t . This tensor is not symmetric. The virtual work of the interior forces takes the first intermediate form W i pu p , w p q ¡ p pu p pp, tq, tq : D p w p ppq dV p0q.For further analyses, it is highly preferable that the weak form of the problem involves symmetric quantities. To address this issue, the second Piola-Kirchhoff stress tensor S p is introduced as S p pu p pp, tqq Jpu p pp, tqq pF p pu p pp, tqqq ¡1 σ x pu x px, tqq pF p pu p pp, tqqq ¡t . pu p pp, tq, tqq t D p w p ppq

					(3.36)
				»
				V p0q	(3.37)
					(3.38)
	The weak form of the interior forces takes the second intermediate form
	W i pu p , w p q ¡	»	V p0q	S p pu p pp, tq, tq :

τ ¡ pF p © dV p0q.

  results in W i pu p , w p q ¡ » V p0q rH p ppq : E p pu p pp, tq, tqs : pF p pu p pp, tq, tqq t D p w p ppq If one develops explicitly H p : E p in terms of derivatives of u p , it leads to H p ppq : E p pu p pp, tq, tq H p ppq : D p u p pp, tq u p pp, tq t D p u p pp, tq One can then separate the linear and non-linear parts in the virtual work of the interior forces by sorting the terms thanks to the properties of the trace operator W i pu p , w p q ¡ p w p ppq : H p ppq : D p u p pp, tq dV p0q p w p ppq : pD p u p pp, tq S p pu p pp, tq, tq dV p0q D p w p ppq : H p ppq : u p pp, tq D p u p pp, tq t © dV p0q. The first term is clearly linear with respect to u p and corresponds to the classical linear elastic interior forces. The two remaining ones are clearly non-linear with respect to u p but all terms are linear with respect to w p . Equations (3.42) to (3.44) are valid for elastic materials only. To deal with damping that requires the use of visco-elastic materials, one can invoke the so-called elastic/viscoelastic equivalence principle [Salençon 1983]. The constitutive equation associated with linear visco-elasticity is quite similar except that (a) tensor H p is replaced by a complex tensor G p and (b) this equation is written in the Laplace domain (variable s describes the "time" dimension). The constitutive equation becomes S p pu p pp, sq, sq G p pp, sq : E p pu p pp, sq, sq.

		»
	¡ ¡	» »	V p0q V p0q V p0q	2 1		(3.44)
	Remark: visco-elasticity
					%	dV p0q. (3.43)
				1 2	H p ppq :

¡ D p D D ¡ D p

  D p δu p pp, tq H p pp, tq : E p pu 0p pp, tq, tq ¨$ : D p w p ppq dV p0q δu p pp, tq t F p pu 0p pp, tq, tq w p ppq t F p pu 0p pp, tq, tq Nppqs t : prH p ppqs : rE p pq 0 qsq : rD p Nppqs dV p0q pq 0 qs t : rH p ppqs : rF p pq 0 qs

		» V p0q	H p pp, tq :	(3.54)
	it is associated with the tangent stiffness matrix, which depends non-linearly on the current discrete displacement tq 0 u, it is symmetric, definite and positive
	rK N L pq 0 qs	» rD p » V p0q V p0q rD p Nppqs t :

¡ D p © : ¡ D p © dV p0q. ¡ rF p © : rD p Nppqs dV p0q.

  ρ p ppq rNppqs t rA yx,¡xy s rNppqs dV p0q

	Ω rD Ω g s, rK c s ¡Ω 2 » V p0q ¡Ω 2 rK Ω c s, ρ p ppq rNppqs t rI xx,yy s rNppqs dV p0q rK a s W Ω » ρ p ppq rNppqs t rA yx,¡xy s rNppqs dV p0q V p0q W Ω rK Ω a s, tf ei u » ρ p ppq rNppqs t r¡Ω 2 I xx,yy W Ω A yx,¡xy s tpu dV p0q V p0q ¡Ω 2 tf Ω ei u W Ω tf W Ω ei u.

  allows to write tq st u ¡Ω 2 rK e s ¡1 tf Ω provided that rK e s is invertible. Neglecting rK c s is here mandatory because it depends explicitly on the rotation speed, which would make the dependence of tq st u not as explicit. This result propagates to any quantity that involves tq st u. In particular, the sum of the centrifugal softening matrix rK c s and the non-linear tangent stiffness matrix rK N L pq st qs becomes rKpΩqs rK c s rK N L pq st qs rK Ω,0 s Ω 2 rK Ω,1 s Ω 4 rK Ω,2 s, : rD p Nppqs t © ¡ rD p Nppqs : tq Ω st u : rD p Nppqs t q : rH p ppqs :

	Ω 2 tq Ω st u,	ei u	(3.66)
			(3.67)
	with
	rK Ω,0 s rK e s ,
	rK Ω,1 s rK Ω c s 4 » V p0q rK Ω,2 s » rD p Nppqs t : rD p Nppqs t : ¡ rH p ppqs : ¡ rH p ppqs : ¡¡ ¡ tq Ω rD p Nppqs : tq Ω st u st u t ©©© ©© : rD p Nppqs dV p0q, V p0q » V p0q rD p Nppqs t : ¡ ptq Ω st u t ¡ rD p Nppqs : tq Ω st u	: rD p Nppqs dV p0q

©©

: rD p Nppqs dV p0q.

  ii.) to supplement other reduction techniques:. The set of vectors tTpΩ j qu, j v1, N Ω w can be • the full set of vectors rT Ω s can be the input of a subsequent reduction process.The orthogonalization procedures are those of the considered reduction technique. Section 3.3.3 presents an example based on the CICR technique introduced in Sec.3.1.2.

1 allow to build a Campbell diagram for each individual Fourier harmonic. It is illustrated in Sec. 3.3.2 in the case of the HEBD model.

(

Table 3 .

 3

	1: Computation time with 57 operating points [BB]
		MMCCM	CCM0	Full
	Operation	Time (s)	Time (s)	Time (s)
	Pre-processing
	Assembly of polynomial	207	207	207
	Multi-Model CCM	599	198	-
	At each iteration
	Resolution	14	4	125
	Total	1 604	633	7 332

  Gradient relative to the coordinates of ¦ x. , .y ¦ (Semi-)inner product relative to operator/matrix ¦ ||.|| ¦ (Semi-)norm relative to operator/matrix ¦ Discrete load tqu Discrete displacement tp qu Discrete harmonic DOF vector tvu Any discrete quantity Matrices rAs Any matrix rcs Observation matrix rDs Coupling/damping matrix res Fourier matrix rEs Unitary Fourier matrix rIs Identity matrix rKs Stiffness matrix rMs Mass matrix rZs Dynamic stiffness matrix rΩs Skew-symmetric rotation speed matrix rΦs Matrix of eigenvectors rTs Set of vectors

	Kronecker product

Vectors

tfu

4 Fourier matrix for real-valued sampled functions

  B.11) From Eqn. (B.10), p y δ is real if δ 0 or δ N s {2 (if N s is even), else p y δ and p y Ns¡δ come in a complex conjugate pair since ds, cosps δ αq cosps pN s ¡ δq αq, sinps δ αq ¡ sinps pN s ¡ δq αq, If one pairs the complex conjugate terms in Eqn. (B.11), it becomes y s pRepp y δ q cos ps δ αq ¡ Impp y δ q sin ps δ αqq . Recall that N s {2 exists only if N s is even. | 2 N s |p y 0 | 2 2 N s |p y δ | 2 N s |p y Ns 2 | 2 . 2 cos pδ αq ¤ ¤ ¤ 2 cos ps δ αq ¤ ¤ ¤ 2 cos ppN s ¡ 1q δ αq 0 ¡2 sin pδ αq ¤ ¤ ¤ ¡2 sin ps δ αq ¤ ¤ ¤ ¡2 sin ppN s ¡ 1q δ αq

					Ns 2 δ0
					2 (B.12)
	Theorem B.1.3. Parseval's theorem
	For a real-valued sampled function, Parseval's theorem becomes
				Ns¡1 ş0	2 ¡1 |y s Ns	(B.13)
	B.1.Another convenient representation involves the Fourier matrix
	res t	" " " " " " " " !	1 . . . . . . 1	1 . . . . . . ¡1	¤ ¤ ¤ ¤ ¤ ¤	1 . . . . . . p¡1q s	¤ ¤ ¤ ¤ ¤ ¤	1 . . . . . . p¡1q Ns¡1

δ1 2

Remerciements

such that y is square-integrable on R (in the sense of the distributions) and let p y be its Fourier transform. Then

The physical meaning of Parseval's theorem is that the energy conveyed by the signal can be computed equivalently either from the signal itself or from its Fourier transform.

|p ypχq| 2 is the power spectral density that has to be integrated over a given bandwith to lead to energy.

B.1.2 Sampled functions

In numerical physics as well as in signal processing, signals are sampled, i.e. any continuous function y is described by the series of N s values y s ypx s q where x s , s v0, N s ¡ 1w, is a series of values of x. When the samples are evenly distributed along x, there exists the relation ds v0, N s ¡ 1w,

x s x 0 s ∆x. One assumes that x 0 0. f s 1{∆x is called the sampling rate. As a result, the Fourier transform of y is also discrete, known at a series of values of χ such that p y δ p ypχ δ q and the values in the series χ δ , δ v0, N s ¡ 1w are also evenly distributed, such that χ δ χ 0 δ ∆χ. with the convention χ 0 0. The sampling frequency is retrived as f s N s ∆χ.

The relation between y n and p y δ defines the direct and inverse discrete fourier transforms. (B.17)

The Fourier matrix res is invertible. Therefore, the operator whose representation is the latter matrix is a bijection. re ¡1 s is explicitely

The Fourier matrix teu is orthognal since the sine and cosine functions are orthogonal in the sense that for any couple pδ, δ I q These orthogonality conditions allow to write L p pV q is the space of functions u defined on V such that

Definition C.2.2. Sobolev space

A Sobolev space W k,p pV q contains all the functions that belong to L p pV q whose derivatives of order 0 ¤ m ¤ k (in the sense of the distributions) belong to L p pV q. Remark 7. In the following, H k pV q will stand for W k,2 pV q. Moreover, H 0 pV q L 2 pV q.

Theorem C.

Inclusion

The Sobolev spaces H k pV q are such that dpk, lq, k l ùñ H l pV q H k pV q.

Proof 3. The proof of this theorem is derived from the definition of the spaces H k pV q.

Remark 8. The previous theorem states that if k increases, the functions are more and more regular.

C.2.2 Inner-products and norms in Sobolev spaces

Theorem C.2.2. M -inner product

x. , .y M,V :

is an inner product.

. is a norm.

Remark 9. ||.|| M,V is commonly denoted ||.|| 2 when H 0 pV q is referred to as L 2 pV q.

Remark 10. Thanks to the properties of inclusion of H 1 pV q in H 0 pV q, x. , .y M,V is also an inner product on H 1 pV q and ||.|| M,V is also a norm on H 1 pV q.

Theorem C.2.4. K-semi-inner product

x. , .y K,V :

is a semi-inner product.

List of Figures