Malte Krack
email: malte.krack@ila.uni-stuttgart.de

Loic Salles
email: l.salles@imperial.ac.uk

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

Vibration Prediction of Bladed Disks Coupled by Friction Joints

The present review article addresses the vibration behavior of bladed disks encountered e. g. in aircraft engines as well as industrial gas and steam turbines. The utilization of the dissipative effects of dry friction in mechanical joints is a common means of the passive mitigation of structural vibrations caused by aeroelastic excitation mechanisms. The prediction of the vibration behavior is a scientific challenge due to (a) the strongly nonlinear and non-uniform contact interactions involving local sticking, sliding and liftoff, (b) topological complexity of the coupled structure, and (c) the multi-disciplinary character of the problem associated with the need to account for structural mechanical as well as fluid dynamical effects. The purpose of this article is the overview and discussion the current state of the art of vibration prediction approaches. The modeling approaches in this work embrace the description of the rotating bladed disk, the contact modeling, the consideration of aeroelastic effects, appropriate model reduction techniques and the exploitation of the rotationally periodic nature of the problem. The simulation approaches cover the direct computation of periodic, steady-state externally forced and self-excited vibra-

tions using the high-order harmonic balance method, the formulation of the contact problem in the frequency domain, methods for the solution of the governing algebraic equations and advanced simulation approaches, including the concept of nonlinear modes. In the ongoing quest for improved efficiency, modern turbomachines are driven near their structural mechanical limits. The successful operation of these machines depends largely on the structural mechanical integrity of the rotating components, owing to their comparatively high loading [152]. Bladed disks undergo high mechanical stress during operation. Static stresses are caused by thermal loads, static fluid pressures and rotation-induced centrifugal loads. Mechanical vibrations, caused by additional dynamic loads of different origins, lead to dynamic stresses. Depending on the static stress level, sustained and high dynamic stresses can lead to high cycle fatigue (HCF). This type of fatigue leads to substantial life cycle costs and presents a major safety issue. A primary goal of the design process is to ensure the structural mechanical integrity. Hence, vibrations are a central concern in the design of aircraft engines as well as industrial gas and steam turbines. (a) Synchronous forced response In the case of forced response, the dynamical loads are caused by the rotation of the blades through the circumferentially inhomogeneous pressure field. Pressure inhomogeneities are caused by aerodynamical blade row interactions and non-uniform inflow conditions. An example for blade row interaction are wakes from upstream stator vanes which leads to the so-called nozzle-excitation. Non-uniform inflow conditions are caused by asymmetries in the flow path due to e. g. struts or casing ovality. In steady operation, the resulting inhomogeneous pressure field is essentially time-invariant in the non-rotating frame of reference. From the perspective of the rotating bladed disk, this pressure field takes the form of a wave traveling with rotor speed. This results in dynamic loading with frequencies being integer-multiples of the rotational frequency, hence the term synchronous excitation. Under the condition of resonance; i. e. , if an excitation frequency coincides with a natural frequency of the structure, the forced vibration response can reach particularly high levels. Resonance coincidences are commonly identified using the Campbell diagram, see Fig. 2.

(b) Flutter Flutter refers to the unstable aeroelastic interaction of a vibrating structure with the surrounding fluid flow. In turbomachinery, flutter is mainly caused by the cascade effect; i. e. , the aerodynamical interference among the blades within a blade row. If this interaction is unstable, the blades receive energy from the unsteady flow as a consequence of their vibration, which represents a self-excitation mechanism. This positive feedback leads to continuously increasing vibration level, until nonlinear effects come into play or the structure fails. Form and frequency of the vibration are usually similar to one of the structure's normal modes of vibration. In contrast to forced response, the oscillation frequency is generally not an integer-multiple of the rotational speed [152,101].

Other vibration mechanisms of aeroelastic nature include vortex shedding and rotating instabilities that occur largely under the condition of partial loading. These unsteady mechanisms cause non-synchronous forced or self-excited vibrations. Besides aeroelastic mechanisms, mechanical effects can lead to vibrations of the bladed disk. Common examples are the rubbing between blades and the casing, vibrations due to torsional or lateral dynamics of the rotor shaft, and extreme events such as ingestion of birds or ice (foreign object damage) or structural failure of individual blades (domestic object damage). Different strategies are pursued to avoid and mitigate vibrations of bladed disks. An important approach is to avoid the excitation of resonant vibrations and aeroelastic instabilities in the operating range. This can be achieved by adjusting the system's dynamical characteristics (e. g. blade counts, natural frequencies) with appropriate design measures. If resonance coincidence cannot be completely avoided, one has to make sure that the vibration response remains within tolerable bounds. This can be accomplished by mitigating the excitation level of the associated wave lengths in the relevant frequency range. An important means to avoid flutter is intentional mistuning, i. e. , a deliberate variation of the blade-to-blade properties. If excessive excitation by resonant forcing or flutter still cannot be avoided for all important modes in the relevant operating range, one may have to increase the damping of the system. The total damping D total of the system is composed as,

Means of vibration reduction

D total = D aerodynamical +D material + D joints D mechanical +D rest . (1)
In the case of resonant forced response, an increase of the total damping can lead to a significant reduction of the vibration levels. In the case of flutter, the aerodynamical damping, D aerodynamical is negative, and additional damping can lead to the complete suppression or, if nonlinear effects become important, the stabilization of vibrations in so-called limit cycles. Mechanical damping, D mechanical , can be grouped into material damping, D material , which is usually comparatively small, and the damping due to mechanical joints, D joints . Dissipation can also have non-aerodynamical and non-mechanical, e. g. electromagnetic origin. Such dissipation mechanisms are accounted for in D rest . The most important sub-group of joint damping, D joints , is friction damping, which refers to the dissipative effects related to dry frictional local sliding in mechanical joints. Friction damping is certainly the most established damping technology of bladed disks. A major drawback of friction damping is that it comes at the cost of wear effects. It should be remarked that many damping technologies with successful applications in other fields, cannot cope with the harsh environment (high temperatures, high centrifugal stresses, corrosive gases), and the strictly limited design space. Also, the application of active or semi-active vibration control strategies is hampered by the requirement of fail-safe operation. Noteworthy alternatives to friction damping include piezoelectric shunt damping [START_REF] Hohl | Active and semiactive Vibration Damping of Turbine Blades with Piezoceramics[END_REF]178], eddy current damping [START_REF] Laborenz | Eddy Current Damping: A Concept Study for Steam Turbine Blading[END_REF][START_REF] Laborenz | Eddy Current Damper for Turbine Blading: Electromagnetic Finite Element Analysis and Measurement Results[END_REF], viscoelastic material damping of coatings [177,[START_REF] Ivancic | Experimental Considerations for Determining the Damping Coefficients of Hard Coatings[END_REF], and impact or particle damping. Friction damping takes place in mechanical joints that are either inherent to bladed disks, such as the ones between blades and roots, or can be introduced additionally, e. g. in the form of underplatform dampers. Some of the most common forms of mechanical joints are illustrated in Fig. 3. Besides damping, additional joints also increase the elastic coupling among adjacent blades. This changes the structure's modal characteristics and is in fact sometimes the primary motivation for the introduction of these joints. By means of an appropriate design of mechanical joints, a considerable miti-gation of vibrations and resulting dynamic stresses can be achieved. This, in turn, leads to an improved structural reliability and decreases fatigue-related costs. Furthermore, this can lead to an increased feasible blade design space and an extended range of tolerable operating conditions, and, thus, contributes to an increased efficiency of the turbomachine.

Scientific complexity of the topic

The design of bladed disks with mechanical joints relies on a profound understanding of the relevant physical phenomena and adequate tools for the assessment of the structural dynamic characteristics. Important characteristics in this context are the vibration level, the resonance frequencies and the mechanical damping. During the design phase, these measures are determined largely by means of vibration prediction, whereas tests are many carried out for validation. The vibration prediction is particularly difficult due to the following aspects.

(a) Nonlinearity The contact interactions in mechanical joints represent strongly nonlinear phenomena. This nonlinearity leads to a strong coupling of different time and length scales. More specifically, the local stick, slip and liftoff phenomena in the mechanical joints, occurring on relatively short time and length scales, have considerable effects on the global vibration behavior of bladed disks, occurring on much longer time and length scales, and vice-versa. Moreover, the dependence on the vibration level needs to be taken into account in the analysis of structural dynamic characteristics, such as resonance frequencies, effective damping and deflection shape. Also, nonlinearity can give rise to phenomena such as co-existence of multiple stable vibration states, and steady-state vibrations that exhibit significant frequency components not present in the excitation spectrum. Suitable simulation methods are often based on iterative, numerical procedures which are comparatively time-consuming.

(b) Model order Turbomachinery bladed disks exhibit blades with generic, three-dimensional profiles and often consist of a large number of components, possibly including additional devices such as friction dampers. Hence, spatial discretization is commonly carried out using finite elements. Moreover, the different components have a number of extended contact interfaces where nonlinear contact interactions may take place. A fine spatial discretization is required to accurately resolve the local contact interactions and the dynamic stress field. This leads to a comparatively high order of the mathematical model.

(c) Multi-disciplinary character : Since the most important vibration mechanisms are of aeroelastic type, both the structural mechanical and the fluid dynamical domain need to be taken into account. These domains can in general not be regarded as independent of each other. For instance, the aeroelastic interaction is essential to understand the physical phenomenon of flutter [152,101]. Besides aerodynamics and structural mechanics, tribology is another scientific field inherently associated with friction damping.

In addition to these aspects, several system parameter are considered uncertain, and may have to be accounted for using probabilistic methods. Finally, many turbomachines undergo a large range of operating conditions. The structural mechanical and aerodynamical properties can vary considerably with the operating condition, which has to be considered during the vibration prediction. The development of suitable vibration prediction methods encompasses a variety of scientific topics. A non-exhaustive list is illustrated in the mind map in Fig. 4.

Motivation and scope of this review article

The development of approaches for the vibration prediction of bladed disks coupled by friction joints was and still is the objective of intensive research efforts. Since the first contributions in the 1970ies [START_REF] Earles | A linearized analysis for frictionally damped systems[END_REF], the degree of detail in the modeling and simulation approaches has gradually increased and reached a considerable level. Considering the sheer number of incremental contributions to the topic, the authors are convinced of the need for an overview. The purpose of this article is the concise presentation and discussion of state-ofthe-art methods for the vibration prediction of bladed disks coupled by friction joints. Instead of an exhaustive presentation of all previously used variants, we present the central notions of general approaches and then focus on established variants. Moreover, it is our intention to provide a thorough presentation and discussion of the typical assumptions underlying the specific modeling and simulation approaches. It should be noted that this goes beyond the scope of most original research articles, although this is essential for a profound understanding of the approaches. Finally, we propose a consistent notation, which is used throughout this article, and we hope that this notation is adopted by the scientific community. In contrast, we do not provide an introduction to the topics of vibrations of bladed disks or friction damping in general. These aspects are extensively covered by the available literature, see e. g. [152]. Useful concepts related to these topics, such as the Campbell diagram, the nodal diameter diagram, and typical characteristics such as damper optimization and performance curves, are therefore not the focus of this work. Previous efforts should be mentioned [126,[START_REF] Firrone | Modelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks[END_REF]135], which also fall into the category of overviews, but aim at presenting rather specific methodologies with interesting numerical examples [126,[START_REF] Firrone | Modelling Friction Contacts in Structural Dynamics and its Application to Turbine Bladed Disks[END_REF], or focus more on the historical development without going into the mathematical details of the approaches [135]. The prospective primary reader group are engineers and researchers in the field of structural mechanics of bladed disks, who might consider this article as a useful introduction to the topic of vibration prediction. They will find this article particularly appealing if they are using or developing according simulation tools. This article also provides an overview for related scientific communities, such as the aerodynamics and aeroelastics communities. Finally, it should be stated that the approaches are presented in a generality that permits the application to the much broader problem class of jointed structures. The article is organized as follows: The modeling is addressed in Section 2. Here, the equations governing the vibrations of a bladed disk, coupled by mechanical joints, and exposed to the fluid flow are derived. Starting from the resulting set of nonlinear ordinary differential equations (ODEs), appropriate simulation methods are addressed in Section 3. In the focus are efficient computation methods for periodic, steady-state vibrations in the dynamic regimes of primary importance, namely synchronous forced response and flutter. An overview of the overall vibration prediction process, and possible directions of future research end this article in Section 4.

Modeling

The bladed disk illustrated in Fig. 5 serves as notional model in this section. The bladed disk rotates at the speed Ω rot and is exposed to a fluid flow. The structure is essentially rotationally periodic; i. e. , the full bladed disk consists of n s almost identical sectors. The sector with index (0) is referred to as the reference sector. Blades and disk can be individual bodies coupled by mechanical joints, or they are integrated in one body as a blade integrated disk (blisk). A sector may also contain one or multiple friction dampers, which are separate bodies. In general, a mechanical joint couples either certain parts within a sector (intra-sector coupling), or Fig. 4: Mind map of the scientific topics relevant for the development of suitable tools for the vibration prediction of bladed disks coupled by friction joints Fig. 5: Schematic illustration of a rotating bladed disk with mechanical joints, possible contact interfaces are indicated in red certain parts of adjacent sectors with each other (intersector coupling). An example for intra-sector coupling is the clamping of the blades in the disk. An example for inter-sector coupling are shroud joints. Depending on the local behavior in the contact interfaces under the considered operating conditions, it is appropriate to describe the mechanical coupling (in specific regions) in either a linear or a nonlinear way. This will be detailed later in Subsection 2.2.

Since the model in Fig. 5 only contains a single bladed disk, it is a priori limited in three respects:

(a) The mechanical contact occurring between the bladed disk and the casing is not considered. Such strongly nonlinear contact interactions between rotating and fixed components are usually temporary in nature, and they require designated modeling and simulation techniques. For a recent overview, the reader is referred to [START_REF] Jacquet-Richardet | Rotor to stator contacts in turbomachines[END_REF]. (b) The mechanical coupling between different stages is excluded. This coupling is usually linear in nature, and it can have considerable effects on the vibration behavior in certain configurations [START_REF] Bladh | Effects of Multistage Coupling and Disk Flexibility on Mistuned Bladed Disk Dynamics[END_REF]. The approaches presented in this article can be modified to take into account bladed disks of multiple interconnected stages [START_REF] Laxalde | Dynamical analysis of multi-stage cyclic structures[END_REF][START_REF] Laxalde | Modelling and analysis of multistage systems of mistuned bladed disks[END_REF]. Hence, the multi-stage coupling does not represent an important conceptual difficulty. We still preferred to exclude the aspect of multi-stage coupling, since it is only rarely considered in the literature and would make the notation of the equations more complicated. (c) The third aspect is related to rotor-dynamic effects.

The bladed disk can generally be excited by the motions of the shaft, and the vibrations of the bladed disk may feedback on the dynamics of the shaft.

To the authors' knowledge, however, this aspect has only rarely been investigated in the context of nonlinearly-coupled bladed disks [START_REF] Gruin | Nonlinear dynamics of a bladed dual-shaft[END_REF], and is therefore not considered in this paper.

This section is divided into four subsections. The modeling of the rotating structure is presented in Subsection 2. Modern turbomachinery blades exhibit three-dimensional profiles, optimized to achieve the desired aerodynamical efficiency. Thus, the spatial discretization of bladed disks is commonly conducted in terms of finite elements.

As mentioned earlier, the structural dynamic characteristics significantly depend on the operating conditions [START_REF]Stationäre Gasturbinen, nachdr. edn. VDI-Buch[END_REF]. It is common practice to regard these conditions as fixed, and thus to assume centrifugal, thermal and static fluid pressure loading as constant. To determine the effect of the static loading, a static finite element analysis (FEA) step is carried out. The static, and in particular the centrifugal loading often results in large deformations of the blades. Mechanical joints that are initially open, may thus close beyond a certain rotational speed. The static FEA thus needs to take into account nonlinear geometrical effects and the nonlinear behavior in the contact interfaces.

Remark The consideration of a fixed operating point is a major simplification. It should be recalled that an important vibration mechanism is the synchronous forced response, where the excitation frequencies are integer multiples of the rotational speed. Fixed operating conditions imply a fixed rotational speed and hence constant excitation frequencies. In a typical forced response analysis, it is, however, crucial to vary the excitation frequency in the neighborhood of specific resonance frequencies. Varying the excitation frequency while regarding a model derived for fixed rotational speed introduces a certain degree of inconsistency. Fortunately, the relevant rotational speed range is typically small, e. g. only a few percent around the nominal speed, so that the assumption of constant mechanical, thermal and aerodynamical characteristics is often justified. In general, the validity of this procedure can be assessed by a point-wise simulation for different rotational speeds. As alternative, the dependence of certain properties on the rotational speed can be explicitly modeled. To this end, an interpolation method is proposed in [153]. This method is, however, so far limited to structures without nonlinear coupling. The motion of the structure is described in terms of nodal displacements measured from the static equilib-rium position in the rotating frame of reference. Small vibrations are commonly assumed. This permits the linear description of mechanical elastic and inertia forces within the structure. This assumption also lays the foundation for the linearization of the aeroelastic forces, which will be addressed later. Consequently, only the contact forces, which are localized at the mechanical joints, are considered nonlinear. This greatly simplifies the structural dynamic analysis. The dynamic deformation of the entire bladed disk, including all blades, disk and, if applicable, dampers and other structural elements is completely described by the nodal displacement vector fe u. The equations governing the system's motion can be expressed as a set of second-order ODEs in fe u,

fe M fe ü(t) + fe D fe u(t) + fe K fe u(t) + fe f c [ fe u] = fe f a ( fe u, fe u, fe ü, t) . (2)
Herein, overdot denotes differentiation with respect to time t. fe M , fe D, fe K are the matrices of acceleration, velocity and displacement proportional structural forces, respectively. Their physical meaning is described in detail below. fe f c gathers the nonlinear contact forces acting on all (intra-and inter-sector) contact interfaces of the bladed disk. The operator notation [•] is used to indicate that, owing to the hysteretic character of friction, fe f c is generally not an explicit function of fe u, fe u at time instant t, but also depends on the time history of these variables. This aspect is further addressed in Subsection 2.2. fe f a denotes the aerodynamical forces and is detailed in Subsection 2.4. The nodal displacement vector, fe u, is organized sectorwise,

(0) fe u T • • • (ns-1) fe u T , where (n)
fe u refers to the nodal displacement vector of sector n. Note that one sector may comprise different individual bodies, typically a disk segment, a blade and possibly a friction damper or other structural elements. The forces fe f c , fe f a are organized in accordance with fe u. The number of degrees of freedom of sector n is n fe,n , such that the total number of degrees of freedom is n fe = ns-1 n=0 n fe,n . The dimension of fe u, fe f c , and fe f a is accordingly n fe × 1, that of the structural matrices is n fe × n fe . If all sectors have identical (aerodynamical, material, geometrical, and contact) properties, the system exhibits ideal rotational periodicity, as opposed to only essential rotational periodicity. For congruent discretization and ordering, the structural matrices take a block circulant form, and the force-deformation relations underlying 1 If, on the other hand, sector-by-sector deviations exist, the system's rotational periodicity is perturbed, which is referred to as mistuning. The mass matrix fe M represents the system's inertia, and it is symmetric and positive-definite. The displacement proportional forces described by fe K comprise linear-elastic forces, as well as stress stiffening and spin softening effects [START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]. Stress stiffening is caused by large static deformations, as a consequence of high static, mainly centrifugal, loading. Spin softening is caused by static deflections in the radial direction, which results in an alteration of the centrifugal forces. As a consequence, fe K generally depends on the rotational speed.

(n) fe f c and (n) fe f a are identical for all sectors n ∈ [0, n s -1].
fe K depends nonlinearly on the static deformation and is obtained by linearization2 . Like the mass matrix, fe K is symmetric, but it is only positive-semidefinite, since rigid body modes may be present. The velocity proportional forces can be split into damping and coriolis forces associated with the symmetric and skew-symmetric part of fe D, respectively. Coriolis forces become relevant if the vibrations have a considerable component perpendicular to the rotation axis, and they depend on the rotational speed. In particular, if blade and disk are relatively compliant, significant radial deflections may occur, which give rise to tangential coriolis forces [171,112]. The effect of coriolis forces on the dynamics of bladed disks is still a matter of research at present, and they have commonly been neglected in the past. We will therefore not further consider coriolis forces in this work. Consequently, fe D is assumed symmetric and reflects only structural damping forces related to dissipative material behavior. In contrast to the fe M and fe K, the damping matrix cannot be derived from conventional FEA. In general, the precise characterization and quantification of material damping forces is still an unresolved problem. In the case of bladed disks, material damping is comparatively weak, but in the absence of other damping mechanisms, it can have a significant influence. Owing to its weakness, material damping forces are usually modeled linearly; i. e. , in terms of a constant damping matrix fe D. The damping matrix is commonly assumed to be symmetric and at least positive semi-definite. It is often defined in terms of logarithmic decrements, obtained e. g. from free-decay experiments. It should be noted that the damping forces do not necessarily have viscous (velocity-proportional) character. The damping term fe D fe u in Eqs. (2) should therefore be regarded as a 'place holder' for any type of linear damping, including linear-hysteretic damping which is very popular in the context of periodic vibrations [START_REF] Ewins | Modal Testing: Theory and Practice[END_REF][START_REF] Gasch | Strukturdynamik. Band 1: Diskrete Systeme, Strukturdynamik[END_REF].

Mistuning

Mistuning refers to the deviation of properties among the sectors within periodic structures. In the case of bladed disks, these deviations can appear as stiffness or density variations within the mechanical structure, geometrical variations of the structure's surface, and variations of the contact properties within the joints. These deviations can be small, for instance in the case of material inhomogeneities, or large, for instance in the case of fatigue cracks or foreign object damages. To a certain extent, these deviations are inevitable due to manufacturing and repair tolerances. During operation, these deviations can be magnified or diminished due to wear and tear. These deviations can affect both the structural mechanical and the aerodynamical behavior of the bladed disk. An important phenomenon associated with mistuning is the localization of the structural dynamic behavior, in particular in the presence of nearresonant forcing. The global vibration is then no longer a (pure) traveling wave, but may typically has a considerable stationary component with a large response in only a few sectors. At these sectors, the vibration response is magnified compared to the tuned case, which is known as mistuning magnification. It is known that the mistuning magnification attains a maximum versus the mistuning level and then becomes comparatively insensitive with respect to additional mistuning. In order to increase the robustness with respect to unintentional and uncertain parameter deviations, it is thus a common strategy to introduce intentional mistuning [START_REF] Castanier | Consideration of the Benefits of Intentional Blade Mistuning for the Forced Response of Turbomachinery Rotors[END_REF][START_REF] Castanier | Investigation of the Combined Effects of Intentional and Random Mistuning on the Forced Response of Bladed Disks[END_REF][START_REF] Lim | Intentional Mistuning Design Space Reduction Based on Vibration Energy Flow in Bladed Disks[END_REF][START_REF] Kenyon | Forced Response of Turbine Engine Bladed Disks and Sensitivity to Harmonic Mistuning[END_REF][START_REF] Kenyon | Maximum Bladed Disk Forced Response From Distortion of a Structural Mode[END_REF]. Mistuning is also known to have a beneficial effect on the aeroelastic stability of blade cascades [START_REF] Kielb | Aeroelastic Characteristics of a Cascade of Mistuned Blades in Subsonic and Supersonic Flows[END_REF]149,164]. The uncertain character of the properties is usually described in terms of its stochastic characteristics. Once these (input) stochastic characteristics are quantified, probabilistic methods, such as Monte Carlo simulations, are utilized in order to determine the stochastic characteristics of the (output) vibration behavior. To estimate the worst or the best configuration with regard to a specific performance measure, optimization techniques are often applied. If the mistuning level is rather small, or the inter-sector coupling is sufficiently strong, the effect of mistuning remains comparatively small [147,116,127,113]. Mistuning is often neglected in the presence of strong coupling by inter-sector contact interfaces. Most of the approaches in this work are kept sufficiently general so that they can be applied to both tuned as well as mistuned systems. We will emphasize the possible simplifications in the tuned case, where appropriate. It should be noted that some of these simplifications are also relevant in the analysis of mistuned systems, since the nominal modes of the tuned system often serve as basis for the mistuned case [START_REF] Lim | A Compact, Generalized Component Mode Mistuning Representation for Modeling Bladed Disk Vibration[END_REF][START_REF] Laxalde | Nonlinear modal analysis of mechanical systems with frictionless contact interfaces[END_REF]. We will, however, not address mistuning-specific methods, since these are extensively covered by the available literature, see for instance the reviews [151, 13].

Traveling wave coordinates

It is generally possible to expand any given physical quantity in terms of its wave components. Due to the cyclically closed topology with only a finite number of n s sectors, the possible wave forms are limited to a discrete set of wave lengths. An exact coordinate transform can be defined between the physical coordinates of the individual sectors and the traveling wave coordinates 3 associated with the different wave lengths. The traveling wave coordinate system represents an apt basis for the description and interpretation of typical vibration phenomena and excitation mechanisms. The mathematical preliminaries and the notion of traveling waves is presented in A. In the following, this coordinate transform is applied to the bladed disk model. For the transformation of the structural matrices A ∈ { fe M , fe D, fe K}, it is useful to distinguish between a block circulant component A and a deviation component ∆A,

A = A + ∆A , A ∈ { fe M , fe D, fe K} . (3)
If each sector has identical matrices, ∆A = 0, and matrices A are strictly block circulant. The transformation of the equations of motion (2) in physical coordinates to traveling wave coordinates yields

tw M tw ü(t) + tw D tw u(t) + tw K tw u(t) + tw f c [ tw u] = tw f a ( tw u, tw u, tw ü, t) -tw ∆M tw ü(t) -tw ∆D tw u(t) -tw ∆K tw u(t) . (4)
In accordance with the transform [START_REF] Gasch | Strukturdynamik. Band 1: Diskrete Systeme, Strukturdynamik[END_REF], the forces are defined as

tw f = tw f T 0 • • • tw f T ns-1 T with tw f k = w H k+1 ⊗ I n fe,s f
, where f denotes either the contact or the aeroleastic forces, f ∈ { fe f c , fe f a }. The application of the coordinate transform [START_REF] Gasch | Strukturdynamik. Band 1: Diskrete Systeme, Strukturdynamik[END_REF] to the block circulant matrices A gives rise to block diagonal matrices tw A = bdiag{ tw A k } on the left hand side of Eq. ( 4), whereas the matrices tw ∆A = W H ns ⊗ I n fe,s ∆A W ns ⊗ I n fe,s are generally fully occupied, where A denotes any of the structural matrices, A ∈ { fe M , fe D, fe K}. The 3 sometimes also referred to as cyclic or nodal diameter coordinates matrices tw M k , tw D k , tw K k can be obtained from the structural matrices of the reference sector by applying appropriate boundary conditions of the form u

(0) l = u (0)
r e iθ k . Since the structural matrices in physical coordinates are assumed symmetric, the structural matrices in traveling wave coordinates are hermitian; i. e. , tw A k = tw A H k . Moreover, since the matrices A are realvalued, the traveling wave structural matrices come in complex conjugate pairs tw A k = tw A j * for corresponding inter blade phase angles (IBPAs) θ k = 2π-θ j , which can be useful for their calculation.

Remark It is important to note that Eq. ( 4) represents a complete description of the full structure, and is thus absolutely equivalent to Eq. ( 2). Only the equivalent discretization of each sector was assumed, which is a prerequisite for the application of the transform defined in Eq. [START_REF] Gasch | Strukturdynamik. Band 1: Diskrete Systeme, Strukturdynamik[END_REF]. So far, no assumption was made regarding the symmetry of the response fe u(t), the aeroelastic forces fe f a or the (typically nonlinear) forces at the contact interfaces fe f c .

Traveling wave coordinates are particularly useful to describe the common aeroelastic excitation mechanisms. As mentioned earlier, in the case of synchronous forced response, the excitation takes the form of a traveling wave from the perspective of the rotating bladed disk. In the case of flutter, usually only the traveling wave forms associated with specific IBPAs are aeroelastically instable. In these cases, the excitation acts on only a small subset of traveling wave components k, whereas tw f a,k = 0 for most wave numbers. The subproblems related to individual wave numbers k are decoupled if (a) the structural properties are identical for each sector ( fe ∆M = fe ∆D = fe ∆K = 0), and if also (b) fe f a and fe f c are linear in fe u and fe u and take the same functional form for each sector. In this case, the nature of the excitation propagates to the response, and hence only the corresponding traveling wave components will exhibit a non-trivial response. This leads to a considerable reduction of the problem and permits to decrease the vibration prediction effort compared to the coupled (even though sparse) n s n fe,s -dimensional problem in Eq. ( 2), especially if n s is large. The subproblems in Eq. (4) related to individual wave numbers k are generally coupled in the mistuned and/or the nonlinear case. The perturbation of the rotational periodicity by mistuning can impede the development of traveling wave vibrations and cause localized vibrations, as discussed in the previous subsubsection. Often, the sector-to-sector deviations are small; i. e. , |∆A| A , such that Eq. ( 4) represents an apt basis for the ap-plication of perturbation techniques. Nonlinearity may also lead to non-trivial response in not directly driven wave components, as it introduces coupling among the different wave components. In particular, higher harmonics (in time and in space) are excited by the action of the nonlinear force. In certain situations, nonlinearity alone may change the nature of response from traveling waves to standing waves, associated with the localization of vibration energy [167,[START_REF] King | A Very Complicated Structure of Resonances in a Nonlinear System with Cyclic Symmetry: Nonlinear Forced Localization[END_REF][START_REF] Georgiades | Nonlinear Modal Analysis and Energy Localization in a Bladed Disk Assembly[END_REF][START_REF] Grolet | Free and Forced Vibration Analysis of a Nonlinear System With Cyclic Symmetry: Application to a Simplified Model[END_REF]. But this represents a strongly nonlinear behavior related to breaking of the underlying symmetries. In the tuned and weakly nonlinear case, the traveling wave character of the excitation is typically propagated to the response. This greatly simplifies the evaluation of the nonlinear force term and permits the reduction of the nonlinear problem to a reference sector with appropriate boundary conditions, as will be shown in Subsection 3.3. Finally, it should be remarked that the formulation in Eq. ( 4) can only be rarely found in the literature. This is because the governing equations of motion are often further transformed to the frequency domain. The time domain formulation given in Eq. ( 4) is considered particularly useful for the investigation of transient vibration problems, such as resonance passages, in the linear or weakly nonlinear regime, for rotationally periodic structures or those with only small mistuning effects.

Most of the approaches presented throughout this work can be applied both in physical and traveling wave coordinates. In order to keep the notation general, we do not explicitly denote the indices fe and tw in the following, such that e. g. u can stand for either fe u or tw u. We highlight the specific features related to either coordinate system, where the distinction is relevant.

Contact modeling

For bladed disks, contact modeling involves describing the essential interactions between dry, rough, usually metallic surfaces of elastic bodies. The bodies are commonly considered as closed systems in the sense that no material transport occurs, as in the case of wear or diffusion. Furthermore, it is usually assumed that the contact interactions do not affect the temperature field such that only the mechanical interactions are taken into account. These mechanical interactions can generally be split geometrically into the normal contact, related to the normal direction of the surfaces in contact, and the tangential contact, related to the interactions perpendicular to the normal direction. For the normal contact, unilateral interactions are commonly considered, and frictional interactions for the tangential contact. Even in the comparatively narrow field of friction damping, it has to be stated that there is not yet a scientific consensus regarding the accurate modeling of contact interactions, and this remains to be a lively field of research. General overviews on the topic of contact modeling can be found e. g. in [174,[START_REF] Johnson | Contact Mechanics[END_REF]172]. In the following, we discuss the different aspects of contact modeling in the light of friction damping of bladed disks. These aspects are grouped into splitting static and dynamic contact problems, contact discretization (kinematics, kinetics), contact laws, and solution approaches for the contact problem.

Splitting the contact problem into static and dynamic subproblems

On the one hand, it is essential to account for the large deformation of the blades due to static forces, which has an important influence on the contact interactions. On the other hand, it is commonly assumed that the vibrations only lead to small deformations. From a mathematical perspective, it is therefore useful to divide the contact problem into a static and a dynamic one. From the solution of the static problem, three results are obtained that are important for the dynamic contact problem: (a) the deformed shape, (b) the active contact regions, and (c) the contact situation in these regions in terms of pressure and clearance distributions. The active contact regions in (b) are defined in such a way that these cover at least those regions that may change their contact behavior due to vibrations, that is, they might liftoff, come into contact, or undergo frictional sliding during contact. More specifically, areas with a large gap do not have to be further considered in the dynamic contact modeling. The same applies to contact areas with excessive normal pressure that will remain in sticking conditions in the presence of vibrations. For the dynamic contact problem, small deformations relative to the static equilibrium position are assumed. Hence, contact detection is usually not a relevant task in the dynamic problem. When the static equilibrium changes, e. g. if a different operating point and rotational speed is considered, the deformation of the contact interface and the initial contact situation generally needs to be updated.

Remark It should be emphasized, that the dynamic variation of the effective contact area is allowed for, which occurs if initially closed (open) contact regions undergo partial liftoff (contact) during vibration. Also, it should be remarked that the static equilibrium position of the nodes is generally not identical to their mean position in the dynamic problem, since the mean gap distribution depends on the vibration behavior. In this sense, the splitting into static and dynamic subproblems is fully consistent with the small vibrations assumption and does not introduce any further approximation.

Contact discretization

It is state of the art to discretize the contact problem in accordance with the finite element description of the underlying continua. The most accurate, but, at the same time, most involved discretization is a Mortarlike 4 formulation which accounts for the consistent coupling of subdomains with generally non-matching meshes [174,172]. The simplest case, on the other hand, is a lumped formulation where the displacement of the interacting surfaces is constrained to certain rigid body motions, controlled by one reference node per surface. A phenomenological generalized force-deformation relation can then be introduced between the reference nodes. In between these two extreme cases, a variety of approaches with intermediate levels of discretization was used in the field friction damping of bladed disks. It should be emphasized that a finer discretization has certainly the potential of a higher degree of predictability. However, a finer discretization alone does not necessarily lead to a more accurate description of the real behavior. In order to assess the discretization quality, it is generally advisable to investigate the convergence behavior of the sought results with respect to an increasing level of discretization. The discretization should be as coarse as possible to avoid spurious computational burden and as fine as necessary to achieve sufficient convergence among the quantities of interest. In the case of friction damping, the measures of interest are typically resonance frequencies and global vibration levels. Often, a relatively coarse discretization is sufficient to determine these global vibrational quantities. In contrast, a much finer local discretization is required to accurately resolve the contact stress field. In general, it should be remarked that such convergence studies are only rarely found in the literature, and, hence, the appropriate contact discretization remains an open question in the modeling of friction-damped systems. Finally, a posteriori error estimates should be mentioned as a potential alternative to conventional convergence studies, see e. g. [START_REF] Coorevits | A posteriori error control of finite element approximations for Coulomb's frictional contact[END_REF]. However, such approaches have not been applied to blade vibration problems to the authors' knowledge.

Fig. 6: Contact discretization

We would like to keep the following mathematical descriptions as general as possible. To this end, we use the notion of contact points, where the contact laws are locally evaluated. These points could be the reference nodes in the case of a lumped formulation, or integration points of contact segments in the case of a Mortar-like formulation. We further assume fully threedimensional contact, since special cases such as pure frictional contact with constant normal load can be easily derived from this. The contact gap, c g k , for the contact pair k refers to the (three-dimensional) relative displacement between the two contact points P k1 and P k2 (belonging to surface 1 and 2, respectively), see Fig. 6. The gap can be split into a scalar normal contact gap g n,k , and a two-dimensional tangential gap vector g t,k ,

c g k = g n,k g t,k , c g =    c g 1 . . . c g nc    . (5)
The normal gap is defined in such a way that a positive value indicates separated contact conditions. The relative displacements are assembled for all n c contact gaps in the vector c g of dimensions 3n c ×1, as indicated in Eq. [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling[END_REF]. Of course, the contact points may belong to disjunct contact interfaces. In accordance with the small vibrations assumption, there is a linear relationship between the contact gaps c g and the coordinates u5 , which can be expressed as,

c g = B T u , (6)
in terms of the constant coupling matrix B of dimension n d × 3n c . The coupling matrix has the form

B = b 1 • • • b nc with the columns b k = b k,1 -b k,2
can be split into portions associated with contact point P k1 and P k2 , respectively. In the case of a node-to-node contact formulation (with conforming meshes), the matrix B is a signed boolean matrix [START_REF] De Klerk | General Framework for Dynamic Substructuring: History, Review and Classification of Techniques[END_REF], provided that the local coordinates systems are aligned with the contact coordinate systems,. In general, B accounts for the possible mapping between non-conforming finite element meshes and other coordinate transformations 6 .

In the case of traveling wave coordinates, the phase lag boundary conditions must be considered. Suppose that contact gap k refers to an inter-sector coupling. In this case, the phase lag boundary condition for wave number j is incorporated into b k by multiplying either b k,1 or b k,2 , whichever is attached to the lower-numbered sector, with e iθj . To obtain the contact forces, the local contact pressure C p (x) = p n p T t T needs to be integrated over the (union of) contact surfaces C,

f c = C C b (x) C p (x) dA ≈ nc k=1 C b (x k ) C p (x k ) ∆A k λ k = Bλ . (7)
Herein, C b (x) denotes the union of shape functions restricted to C, which relates the nodal coordinates assembled in u with the relative displacement field at the contact interface as C g (x) = C b T (x) u. As indicated in Eq. ( 7), the continuous integral is approximated by a weighted sum over a finite set of n c contact points, where x k and ∆A k refer to the location and the area associated with contact point k, respectively, see Fig. 6. The contact laws presented in the following define a relation between C p and C g (and their time derivatives). Note that impact-type contact behavior could be accounted for equivalently to Eq. [START_REF] Bonhage | Transient amplitude amplification of mistuned blisks[END_REF], where the forces had to be replaced by impulses.

Remark In order to ensure consistency with finite element modeling, the contact laws generally have to be evaluated on the pressure level instead of the force level. Otherwise, the different weights ∆A k are neglected, which leads to inconsistencies in the case of non-regular grids or at the boundary of regular grids.

Note that the coupling matrix B also relates the local with the global kinetics. In fact, it might be easier to determine B from Eq. ( 7) if a finite element approach is pursued. In the case of a lumped formulation, the integral and shape functions in Eq. ( 7) should be considered notional, since B can usually be directly formulated in this case. In practice, it is useful to break down Eq. ( 7) into the different contact interfaces, and to formulate the related sub-matrices of B only for the nodes directly involved in the respective interface discretization.

Contact laws

Contact laws define relationships between the local contact kinetics (pressure or force) and the local kinematics (gap and/or gap velocity). Throughout this work, we use pressure-based formulations, while force-based formulations can be obtained by spatial integration using Eq. [START_REF] Bonhage | Transient amplitude amplification of mistuned blisks[END_REF]. Besides pressure/force laws, impact laws are commonly used, particularly in the field of rigid body dynamics. It is still an open question if it is reasonable to apply impact laws to elastic bodies described in terms of finite elements. To the authors' knowledge, impact laws have never been used in the context of bladed disks coupled by friction joints. Hence, we exclude impact laws from the further discussion. In general, the contact laws should be as detailed as necessary, to account for the expected interactions in the dynamic regime of interest. On the other hand, they should be as simple as possible; in particular, it is desirable that they only involve a small number of empirical parameters to be obtained from measurements or experience without loosing the main tribological characteristics of the surfaces in interaction.

Normal contact For the normal contact, a unilateral law is commonly used which restricts the interpenetration of the interacting surfaces and avoids the transmission of tensile forces (absence of adhesive effects). Due to vibrations, the normal contact generally undergoes oscillating normal load, including the possibility of vanishing normal load in the case of separation. The dynamic contact interactions depend on the static situation; if the contact is initially open, the contact may close due to vibrations. Similarly, if the contact is initially closed, liftoff is possible if the vibrations are large enough. The normal load determines the limit friction force and thus has a considerable influence on the tangential contact interactions. Unilateral contact represents a strong nonlinearity and may significantly affect the vibration behavior, even in the frictionless case. Different unilateral normal contact laws are depicted in Fig. 7a. Apart from the non-regular rigid formulation, the linear-elastic unilateral law is widely used. We discuss rigid vs. compliant contact laws at the end of this subsubsection. Moreover, nonlinear relations may be employed to account for the effect of varying effective contact area, e. g. due to the deformation of the local surface roughness asperities [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF]172,174]; see [146,159] for applications to bladed disks. The normal contact interaction is typically considered conservative, although few approaches can be found that model the dissipative character of normal impacts in terms of a coefficient of restitution or viscous damping elements.

The rigid and linear-elastic unilateral laws are particularly popular in the field of contact modeling of bladed disks. Also, the treatment of these two laws in the context of the harmonic balance method is also be addressed in Subsubsection 3.3.5. Therefore, we recap their mathematical formulations. For the rigid formulation, the contact law assumes the form of a linear complementarity inequality,

p n = 0 ∧ g n -g n,0 ≥ 0 separation p n ≤ 0 ∧ g n -g n,0 = 0 contact . (8)
For the elastic formulation, the contact law is regular and can be written in explicit form,

p n = 0 g n -g n,0 ≥ 0 separation k n (g n -g n,0 ) g n -g n,0 < 0 contact . (9)
Herein, k n is the normal stiffness per area value.

Tangential contact The interface behavior can generally be modeled using dedicated constitutive approaches. It has been proven that discontinuous behavior like friction [107, 108] could be modeled by means of the free energy and a specific pseudo-potential including thermal effect and wear phenomena. This procedure allows to formalize rigorously the definition of the interface behavior in accordance with the thermodynamic laws. Strmberg [154,155,156] proposed defining interface laws by taking into account the thermodynamic principles as well as all possible wear mechanisms. We will not repeat the theoretical developments at this point, but instead provide a clear framework for the interface behavior using specific rheological models. In this model, the different regimes of sticking and sliding behavior between the dry surfaces are commonly distinguished in the tangential contact plane. Friction has dissipative character and is the primary cause for damping in mechanical joints.

Common friction laws are illustrated in Fig. 7b in terms of their steady-state hysteresis for harmonic, onedimensional motion g t and constant normal pressure p n . It should be noted that p n is generally not constant, as explained in the previous paragraph. For the general case of spatial friction, the two-dimensional character of the relative motion needs to be taken into account.

It should be noted that it is common practice to account for each direction in a decoupled manner by treating them as independent, one-dimensional friction laws, see e. g. [124,150,180]. This can significantly decrease the computational burden in the dynamic analysis [START_REF] Cigeroglu | A Microslip Friction Model with Normal Load Variation Induced by Normal Motion[END_REF]. However, substantial quantitative differences between the coupled and the uncoupled approach are generally possible [125]. A qualitative difference is the possibility of a periodic sliding state which can only exist in the case of spatial friction.

The rigid and elastic Coulomb laws distinguish strictly between sliding and sticking contact. As opposed to this, they do not take into account microslip implicitly. The microslip regime is of particular interest, since it provides a comparatively high dissipation, while retaining a relatively strong elastic coupling between the interfaces. Mircroslip effects are particularly relevant for large rough interfaces of compliant bodies subject to high normal loads [START_REF] Gu | Advanced Modelling of Frictional Contact in Three-Dimensional Motion When Analysing the Forced Response of a Shrouded Blade[END_REF]. When the contact area is discretized and either the rigid or elastic Coulomb law is applied locally, the effect of microslip can be described. A comparatively fine spatial discretization is required to accurately capture the microslip behavior in this case. Instead, or in addition to the fine discretization, microslip effects can be taken into account implicitly using an appropriate hysteretic law, as indicated in Fig. 7b.

Common examples are the Dahl [START_REF] Dahl | Solid Friction Damping of Mechanical Vibrations[END_REF], the LuGre [173], the Bouc-Wen [169], the Iwan [144], the Valanis [START_REF] Iranzad | Identification of nonlinear bolted lap joint models[END_REF] and the Preisach [168] models. It can be generally stated that the microslip laws are frequently applied to contact problems with time-constant normal load. In the case of variable normal load, and in particular in the case of liftoff, numerical difficulties typically arise using these models [125]. In the context of microslip, multi-scale approaches for the contact formulation between rough surfaces should also be mentioned [176], which have, however, not been applied to blade vibration problems to the authors' knowledge. Analogous to the normal contact case, the rigid and elastic Coulomb laws are particularly popular in the field of contact modeling of bladed disks, and we therefore provide their mathematical formulation at this point.

For the rigid formulation, the contact law assumes the form of a nonlinear complementarity inequality, Note that p n ≤ 0, and, thus, p t points into the direction of -ġt during sliding. In the general field of tribology, distinct values of the friction coefficient µ are often used for sticking and sliding, and a characteristic dependence of µ on the relative sliding velocity is frequently used. In the field of blade vibrations, however, a constant value for µ is typically assumed, as is done throughout this work.

ġt = 0 ∧ p t < µ |p n | sticking ġt = 0 ∧ p t = µ |p n | ġt ġt sliding . (10)
In contrast to the set-valued law in Eq. ( 10), the elastic Coulomb law is regular. The hysteretic character of the tangential pressure p t is governed by the differential equation,

dp t = k t dg t p t + k t dg t ≤ µ |p n | sticking 0 p t + k t dg t > µ |p n | sliding . (11)
Herein, k t is a scalar tangential stiffness per area value. This most common case implies isotropic tangential stiffness. To account for anisotropic friction, the scalar k t is replaced by a two-dimensional matrix.

Rigid vs. compliant contact laws

The question whether a rigid or a compliant model better describes the real physical behavior is an ill-posed one. The appropriate choice of the contact model is always strictly related to the appropriate choice of the model of the underlying bodies. Just as it is impossible to tell where surface roughness asperities end and where the underlying body begins, the integral model of both contact and underlying bodies should be regarded as a unit and these two sub-models must be well inter-coordinated.

Typically, the finite element model only accounts for the macroscopic geometry of the contact interface, and the surface roughness is not explicitly captured. An inherent artificial stiffness is introduced by finite discretization. Depending on the resolved length scales, it might be appropriate to neglect this artificial stiffness and to use a rigid contact law. As alternative, one can take the neglected compliance into account in a 'smeared' way by using a compliant contact model. To this end, the normal and tangential stiffness can be utilized to account for the neglected compliance. This represents a physical motivation for using compliant contact models. It should be emphasized that in this case, the parameters k n , k t are related to the discretization, and, therefore, we regard it as inappropriate to refer to them as 'contact stiffness'.

Remark In fact, the idea of the lumped discretization is to deliberately neglect the structural compliance of the region close to the contact interface by considering it as rigid in the model of the underlying bodies, and to account for its compliance instead in the contact model [START_REF] Claeys | Experiments and numerical simulations of nonlinear vibration responses of an assembly with friction joints -Application on a test structure named 'Harmony[END_REF]. In this case, k n , k t mainly represent the structural stiffness. While this approach introduces notable inaccuracies, the resulting model has comparatively good mathematical properties (low to moderate numerical stiffness).

From this discussion, it can be ascertained that the question whether to use a rigid or a compliant contact model boils down to the question if and, if so, where to account for the residual compliance. This question, in turn, is mainly of mathematical nature, as it influences the solution method for the contact problem, as addressed in the following subsubsection. In any case, the modeling approach must be consistent with the experimental parameter identification method.

Solution approaches for the contact problem

The mathematical properties of the contact laws are essential for the choice of suitable solution approaches.

In the case of the elastic formulation, the contact law is regular, giving rise to ODEs, which can generally be solved using standard methods. An important mathematical property of these equations is their numerical stiffness. The ODEs are stiff, if the contact model is stiff (e. g. large values of k n , k t ) compared to the characteristic stiffness of the underlying structure. This can result in ill-conditioning of the problem and the need for a comparatively fine time discretization.

Table 1: Solution approaches and their properties for rigid and compliant contact laws contact law regularization non-smooth formulation rigid by penalty method, standard approach resulting in stiff ODEs compliant already regular, re-formulation using stiffness depends on model supplementary nodes

In the case of the rigid formulation, the contact law is set-valued. The motions are not simply governed by the ODEs, but the complementarity inequalities must be satisfied as well. The governing equations can be stated as a constraint optimization problem involving the Karush-Kuhn-Tucker optimality conditions [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF]. One of the methods qualified for the approximate solution of constrained optimization problems is the penalty method. If a linear penalization with according coefficients is used, this approach becomes mathematically fully equivalent to the elastic formulation. It must be stressed, however, that there is a substantial difference in the physical interpretation of these contact formulations. While k t and k n have a physical motivation and a specific finite value in the case of the elastic formulation, they are the penalty coefficients in the case of the rigid formulation, and thus a mathematical parameter. The penalty coefficients should be specified as large as possible to ensure that the constraints are not significantly violated. Hence, the numerical stiffness dilemma is inevitable in this case. For rigid contact laws, it is more appropriate to employ robust numerical methods specifically designed for non-smooth problems, such as the augmented Lagrangian method. The literature on such methods is rich. In this work, we focus on the solution of the contact problem in the case of periodic vibrations using the framework of spectral methods, and the harmonic balance method in particular, see Subsubsection 3.3.5. Consider the case of a compliant contact model with comparatively high stiffness. In order to overcome numerical difficulties, an interesting approach is to reformulate the problem into a non-smooth one, and use appropriate solution methods. This can be achieved by introducing massless supplementary nodes attached to the contact points via the stiffness characteristics of the contact model [START_REF] Charleux | Étude des effets de la friction en pied d'aube sur la dynamique des roues aubagées[END_REF]. Thus, the comparatively high stiffness is moved to the structural model. At the supplementary nodes, rigid contact laws are considered, rendering the problem non-smooth. Hence, both rigid and compliant contact models can, in principle, be solved using either a regularization or a non-smooth formulation. However, not all combinations are recommendable. The possible solution approaches are summarized in Tab. 1. Besides physical and mathematical reasoning, one may characterize the affinity of certain research groups to either regularized or set-valued formulations as a 'fetish'.

Dynamic substructuring

The number of nodal DOFs of finite element models of a bladed disk can be considerable. In the industrial design process, typical numbers are in the order of magnitude of several millions [122]. With today's computational resources, comprehensive nonlinear dynamical investigations are practically intractable without further model reduction. To reduce the computational burden, dynamic substructuring is commonly applied. Dynamic substructuring of nonlinearly-coupled bladed disks involves two steps: (a) preparing the nonlinear interface coupling for the dynamic analysis, and (b) the reduced description of the underlying linear structure. These two steps are described in the following two subsubsection. Dynamic substructuring can be directly applied to the whole bladed disk in physical coordinates. As mentioned in Subsubsection 2.1.3, the use of traveling wave coordinates is more useful in a rotationally periodic setting. In this case, only a reference sector is considered and appropriate boundary conditions are applied for each relevant wave number. The methods presented in this subsection are applicable to both cases and, hence, the notation is kept general. The specific features of cyclic boundary conditions are explicitly stated, where appropriate.

Preparation of the nonlinear coupling

In the dynamic substructuring procedure, special attention is paid to the treatment of coupling interfaces, at which nonlinear contact interactions are taken into account in the later dynamic analysis. It is therefore useful to distinguish between the coupling coordinates, i. e. , the physical coordinates directly involved in the description of the contact interfaces, and the remaining coordinates which describe those motions that do not directly influence the relative interface deformation. It is common practice to retain the coupling coordinates in the reduced description of the structure, whereas the interior dynamics is approximated in terms of a reduced set of generalized coordinates (by the methods described in the next subsubsection). This strategy permits an accurate description of the local kinematics of the interface. Also this makes it easier to take into account the static residual compliance of the interface, i. e. to retain a statically complete description of the interface [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF]. In case (a), only the interior DOFs, which are not associated to any of the interfaces, correspond to the remaining coordinates. In contrast, in case (b), those interface coordinates that describe the absolute motion, are also excluded from the coupling coordinates. The set relations between the different types of coordinates are illustrated in Fig. 8.

Two variants for the definition of coupling coordinates

The most important benefit of variant (b) is that the number of coupling coordinates is reduced by factor of approximately two7 compared to variant (a). This can significantly decrease the computational effort of the dynamic substructuring procedure, and reduce the number of unknowns considered in the nonlinear dynamic analysis. Furthermore, the use relative coordinates simplifies the coupling formulation, since the contact laws are formulated in terms of the contact gaps. Hence, no further coordinate transformation is necessary in the nonlinear dynamic analysis. Finally, when the relative coordinates are retained, it is easier to locally introduce a problem-adjusted approximation order, for instance by considering a larger number of base functions to accurately resolve the local dynamic behavior of the contact interface. It should be noted that variant (a) is widely used in the vibration analysis of bladed disks with contact interfaces, while variant (b) is much less common. In C, we describe how the transformation to relative coordinates is applied.

Reduced description of the underlying linear structure

Consider a substructure, such as the reference sector of a bladed disk or a sub-domain of it, described either in physical or in traveling wave coordinates. The substructure has a number of n d coordinates assembled in the vector u. The general idea of component mode synthesis (CMS) techniques is to approximate u via a set of n r vectors of component modes, assembled as columns in the matrix T , and associated generalized coordinates assembled in the vector r u,

u ≈ T r u . ( 12
)
Substituting the linear transform defined in Eq. ( 12) into the equation of motion generally results in an error term. This error term is then made orthogonal with respect to the base vectors. This procedure gives rise to the projected equations,

r M r ü(t) + r D r u(t) + r K r u(t) + r f c [ r u] = r f a ( r u, r u, r ü, t) , (13)
where the projected matrices are read r A = T H AT for any of the structural matrices A ∈ {M , D, K}, and the projected force vectors read

r f = T H f for f ∈ {f c , f a }.
Thus, the problem is reduced to the subspace spanned by the component modes, and the number of unknown coordinates is reduced to n r . Hence, the nonlinear dynamic analysis can be solved in the reduced space of component modes, where typically n r n d such that considerable computational savings are often accomplished. In general, any type of base vector could be used such as modal deflection shapes (under appropriate boundary conditions) and static deflection shapes. The various CMS techniques mainly differ in the selection of the component modes. The selection of appropriate base vectors has a crucial influence on the approximation quality and convergence. The number of component modes assembled in T should be as large as necessary to capture the vibration behavior in the dynamic regime of interest with sufficient accuracy, and as small as possible to achieve a significant model order reduction. To assess the accuracy of the selected CMS method, it is common to conduct a convergence study with respect to the number of considered component modes. To this end, correlation measures such as the Modal Assurance Criterion can be evaluated to compare approximated mode shapes with a reference. Most CMS techniques define the component modes based on the underlying undamped structural mechanical problem, that is, by neglecting the damping matrix. The literature is rich in reviews of CMS techniques, see e. g. [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF][START_REF] De Klerk | General Framework for Dynamic Substructuring: History, Review and Classification of Techniques[END_REF]. Remark It is important to note that the basis vectors assembled in the matrix T do not depend on the displacement. Hence, the CMS formalism represents a linear transform and projection procedure. This is reasonable as long as the nonlinearities are local [145,[START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF]. In the presence of global nonlinearities, as in the case of large dynamical deflections of the blades, a variable base in conjunction with a nonlinear projection technique seems more appropriate. Some of the most popular methods Probably the most popular CMS methods are the Craig-Bampton (CB) method and the MacNeal-Rubin (MR) method. Both methods permit to maintain a statically complete description of so-called retained DOFs. Moreover, the transformation matrix T is brought into such a form that these DOFs are still directly accessible in the reduced problem, as opposed to the remaining deleted DOFs. To accurately represent the dynamical behavior, the component modes include a set of normal modes. Mathematical details on the computation of T is given in D for both methods. The CB method is particularly known for its numerical stability and great modal convergence. It is known that the MR method leads to comparatively large computational effort and may suffer from numerical difficulties in some situations [START_REF] Gasch | Strukturdynamik. Band 2: Kontinua und ihre Diskretisierung[END_REF]. Batailly et al. [START_REF] Batailly | Study of Component Mode Synthesis Methods in a Rotor-Stator Interaction Case[END_REF] compared the MR with the CB method for a dynamical contact problem and found that both methods lead to similar results for reasonably large numbers of modes.

Definition of the set of retained DOFs It is common practice to retain not only the coupling DOFs, but also other DOFs of particular relevance, such as DOFs where the vibration response is of interest or DOFs where local forcing is applied. Conceptually, this is not necessary since the forcing can be projected into the reduced space, and the response of any physical DOF can be determined a posteriori by means of expansion using Eq. [START_REF] Castanier | Investigation of the Combined Effects of Intentional and Random Mistuning on the Forced Response of Bladed Disks[END_REF]. In fact, it influences the boundary conditions for the considered normal modes in the case of the CB method, if DOFs other than the coupling DOFs are retained. This tends to degrade the convergence behavior of the reduction basis such that a much larger number of normal modes needs to be retained in order to achieve the sufficient accuracy [START_REF] Charleux | Threedimensional Multiharmonic Analysis of Contact and Friction in Dovetail Joints[END_REF]. Hence, it is recommended to retain only the coupling DOFs in conjunction with the CB and the MR method.

Accurate representation of the static interface behavior

The local relative displacement in the contact interface is typically small compared to the global displacement of the structure. An accurate description of the local elasticity is essential for the modeling of contact interactions. This needs to be ensured by the reduction method. Both the CB and the MR method permit the complete description of the static deformation behavior of the contact interface. In contrast, the modal truncation method in its conventional form is not suitable for contact problems, since the first modal deflection shapes barely contain information about the elasticity of contact interfaces. To describe the dynamic compliance of the coupled structure it is therefore useful to augment the dynamic compliance associated to the considered modes by the so-called static residual compliance, i. e. , the static compliance of the neglected modes [122].

Comment on the use of relative coordinates at the interface If the transformation to relative coordinates at the interface is applied, as defined in Eq. ( 51), prior to computing the reduction basis, only the relative coordinates are retained as coupling DOFs. In this case, the 'fixed interface modes' involved in the CB formulation can be identified as the modes for tied contact conditions, while the interface itself can deform with respect to the inertial frame of reference. In contrast, in the case of the classical formulation, the fixed interface modes represent the modes for literally fixed interfaces; i. e. , where the surfaces are fixed to the inertial frame of reference. This does not represent a physically relevant situation. Since tied contact conditions are an important limit case of the nonlinear contact behavior, it is useful to exactly describe the associated modal deflection shapes. For instance, this property can be helpful for the consistency with boundary conditions imposed in computational fluid dynamics simulations or experiments. If the transformation to relative coordinates is not applied, a comparatively large number of modes is necessary to accurately describe the behavior of a tied interface.

On the number of coupling DOFs Both the CB and the MR method have the disadvantage that a comparatively large number of retained DOFs is required at the interface. More specifically, the number of retained DOFs depends on the level of spatial discretization at the interface. While a comparatively coarse discretization may be sufficient to approximate the relative displacement field at the interface, a much finer discretization is required to accurately resolve the contact stress field. Procedures have been developed to reduce the number of coupling DOFs, see e. g. [161,[START_REF] Donders | A Wave-Based Substructuring Approach for Fast Modification Predictions and Industrial Vehicle Optimization[END_REF]179]. However, in the presence of nonlinearity, a reduction of the interface description is often associated with a significant decrease in accuracy [145,[START_REF] Cook | Concepts and Applications of Finite Element Analysis[END_REF], or a rather limited range of validity.

Specific features in the case of traveling wave coordinates For rotationally periodic sturctures, the CMS procedure is applied to the equations of motions (4) formulated in traveling wave coordinates. Since the traveling wave structural matrices generally depend on the associated wave number, the reduction basis T and the projected structural matrices will also exhibit this dependence. Hence, the CMS procedure needs to be applied for each relevant wave number. Which of the possible n s wave numbers are relevant, depends on the excitation properties and the dynamic regime, as will be explained in Subsection 3.3. As explained in Subsubsection 2.1.3, the traveling wave structural matrices come in complex-conjugate pairs (for real-valued matrices in physical coordinates). This property can be utilized in the dynamic substructuring procedure. Hence, the actual computation of the projected structural matrices must be carried out at most for ns 2 + 1 or ns+1 2 wave numbers for even or odd numbers of sectors n s , respectively.

Cascading the dynamic substructuring procedure

In general, dynamic substructuring can be applied in multiple steps. For instance, the bladed disk is divided into the n s sectors, and each sector is further divided into smaller substructures such as blade, disk section and friction dampers. Reduced models for these substructures are then derived and subsequently coupled to obtain a reduced description for the entire system. This multi-step strategy produces a computational overhead for the preliminary reduction steps, while computational effort for the generation of the individual reduced models is decreased [START_REF] Géradin | Mechanical vibrations: theory and application to structural dynamics[END_REF]. In a cascade-like procedure, the transformation to relative coordinates can be carried out at different stages. An interesting variant is to first reduce the nonlinearlycoupled components, while retaining all DOFs at the coupling surface. Starting from the reduced component models, the transformation to relative coordinates can be applied. Finally, a second CMS step is carried out in which only the relative coordinates at the interface are retained. It should be noted that the transformation to relative coordinates generally has a detrimental effects on the sparsity of the structural matrices. For the computation of the component modes, this means that the linear algebraic operations (eigenvalue analysis, solution of linear equations) are computationally more involved. As mentioned earlier, friction dampers are sometimes approximated as rigid bodies. This can be interpreted as a first dynamic substructuring step, where the elastic modes are neglected.

Consideration of aeroelastic effects

As mentioned in the introduction, two of the major causes for vibrations of bladed disks are of aeroelastic nature, namely forced response and flutter. For the vibration prediction, it is therefore essential to take into account the aeroelastic effects introduced by the surrounding fluid flow. In general, the aeroelastic effects on the structure can be split into external forces f ae , and coupling forces f ai ,

f a (u, u, ü, t) = f ae (t) -f ai (u(t), u(t), ü(t)) . ( 14
)
The external forces explicitly depend on time, but not on the structural vibrations u(t), u(t). Consequently, these forces are also present if the vibrations cease to exist. In contrast, the coupling forces describe the mutual interaction between structural and aerodynamical forces, and therefore depend on the structural vibrations, but they are typically assumed not to have an explicit time dependence. These forces vanish in the absence of vibrations. External forces lead to forced vibrations. Coupling forces mainly introduce aerodynamical damping. If the aerodynamic work done per cycle (or equivalently, the damping) is negative, for a particular form of vibration, an aeroelastic instability is present. As a consequence, the coupling forces can lead to self-excited vibrations (flutter). Moreover, the coupling forces introduce elastic coupling among the blades in a blade row, which leads to a shift in the natural frequencies, and alters the modal deflection shapes. However, the latter effect is often neglected considering the typically large mass ratio between the structure and the fluid. Even in the absence of flutter, the aerodynamical damping is of foremost importance for forced response analyses, since the vibration level in the resonant case is largely controlled by this quantity.

External forces

Nature of external forcing The causes for synchronous aerodynamical forces were already described in the introduction. The frequencies of this forcing are integer multiples of the rotational speed Ω rot . Hence, this form of excitation is termed engine order excitation. A characteristic feature of engine order excitations is their symmetry order around the circumference, that is, their so-called engine order m where m ∈ N + positive integer. An example for a distribution of the pressure p around the circumference is illustrated in Fig. 9, which corresponds to a dominant engine order m = 3. Note that the steady fluid pressure field, which could be described by a zeroth engine order, is already taken into account in the static analysis, cf. Subsection 2.1. In the case of nozzle excitation, this number corresponds to the number of guide vanes of the upstream vane cluster, and is thus comparatively large. Lower engine orders are typical for casing ovality and inlet or outlet asymmetries.

For a constant rotational speed Ω rot , an engine order excitation of specific order m introduces a periodic external forcing of the bladed disk with frequency Ω = mΩ rot . Due to rotation, adjacent blades are exposed to the same pressure with a constant time delay. Mathematically, the pressure (n) p(t) at a certain location on the surface of blade n, at time t, is thus related to the pressure (0) p(t) at the corresponding location Fig. 9: Engine order excitation with dominant engine order m = 3 on the reference blade, (n) p(t) = (0) p(t + n∆t), with ∆t = 2π mΩrotns . Hence, the forcing takes the form of a traveling wave from the perspective of the rotating bladed disk 8 . If multiple engine orders are present, the periodic excitation can be expanded in a Fourier series,

fe f ae (t) = m∈M fe F ae,m e imΩrott =             1 e i 2π ns m . . . e i 2π ns mns       ⊗ I n fe,s       (0) fe F ae,m e imΩrott . ( 15
)
Herein, fe F ae,m is the forcing amplitude of the full bladed disk associated with the engine order m. Often only few of the engine order components fe F ae,m exhibit a considerable magnitude and have to be considered in the simulation. The engine orders m associated with significant components are collected in the set M. Owing to the traveling wave character, the forcing of the full bladed disk can be expressed by the forcing (0) fe F ae,m of the reference sector, as indicated in Eq. ( 15). This property can be utilized to simplify the numerical computation of the excitation forces. For the nonlinear vibration prediction, it is useful to define the so-called fundamental engine order m 0 , as the greatest common divisor, m 0 := gcd M, among the engine orders within the set M. The form of the engine order excitation, as defined in Eq. ( 15), is well-suited for the description in the traveling wave coordinate system. To this end, the transform defined in Subsubsection 2.1.3 is applied to Eq. ( 15). It can be ascertained that a specific engine order component, m, results in only a single non-trivial traveling wave component tw f ae,k ,

tw f ae,k (t) = √ n s m∈M k (0) fe F ae,m e imΩrott , (16)
with

M k = {m | m mod n s = k, m ∈ M}.
Numerical computation Computational fluid dynamics tools, based e. g. on the unsteady Reynolds-averaged Navier-Stokes (URANS) equations, are commonly utilized to determine the excitation forces in terms of the engine order components

(0)
fe F ae,m . Typically, a unidirectional approach is followed by prescribing only the rigid body rotation, but not the dynamic deformation of the blades [START_REF] Filsinger | Approach to Unidirectional Coupled CFD-FEM Analysis of Axial Turbocharger Turbine Blades[END_REF][START_REF] Green | Scaling of Turbine Blade Unsteady Pressures for Rapid Forced Response Assessment[END_REF]163]. Depending on the major source of the flow defects, the system boundary for the fluid dynamic computations is defined. Appropriate boundary conditions need to be imposed at the inlet and the outlet of the considered domain. Moreover, the data transfer between rotor and stator must be ensured, e. g. by using the mixing plane or the sliding mesh technique [163]. Finally, the fluid pressure at the blade surface needs to be mapped from the fluid domain to the finite element model. The time-dependent fluid pressure field p(x, t) is then integrated over the blade surface (0) S b of the reference sector to determine the consistent nodal forces. The Fourier transform is applied to determine the different engine order components m ∈ M,

(0) fe F ae,m = 1 2π (2π) (0) S b η S (x)p(x, t)dA e -imΩrott d (Ω rot t) . (17)
Herein, η S (x) is the vector of finite element shape functions associated with the vector of nodal displacements (0) fe u, restricted to the surface normal direction (with S u(x) = η T S (x)u) 9 . The discretization of the surface integral in Eq. ( 17) leads to a linear relationship between the vector of nodal forces and the vector assembling the pressure evaluated at the discrete integration points (e. g. the nodes of the underlying finite element model). The formulation in Eq. ( 17) is limited to the effect of the fluid pressure. Shear forces can be accounted for analogously. The forces can then be projected onto the reduced basis T .

Remark For the actual implementation, it is useful to distinguish between steady and unsteady pressure. The steady fluid pressure does not contribute to the dynamic loading, provided that η S (x) is constant in time, i. e. , the structure's does not undergo dynamic deformation, which is a common simplification in forced response investigations.

Linearized coupling forces

Theoretical concept It is the current state of the art in the vibration prediction of bladed disks coupled by friction joints to describe the aeroelastic coupling in a linear form. Moreover, it is commonly assumed that the mechanisms underlying the aeroelastic coupling are of rotationally periodic nature. Hence, the aeroelastic coupling between different traveling wave forms is neglected, and the coupling is defined separately for each relevant IBPA. The coupling forces are usually expressed in a displacement-proportional form in terms of a complexvalued, aeroelastic transfer matrix 10tw G ai,k ,

tw f ai,k = tw G ai,k (Ω) tw u k . (18)
Herein, Ω is the oscillation frequency. The coupling forces as defined in Eq. ( 18) can then be substituted into Eq. ( 4) or, more specifically, into its projected form with respect to the component modes.

Often the transfer matrix is defined with respect to a given set of in-vacuum normal modes of vibrations for appropriately linearized contact conditions. The transfer matrix is then typically diagonal dominant, which means that the modal deflection shapes are not significantly affected by the surrounding fluid flow. In this case, the imaginary part of a specific diagonal element of tw G ai,k corresponds to the aerodynamical damping of the associated mode, while the real part is related to the shift in the natural frequency caused by the aeroelastic coupling. The classical (linear) flutter analysis, the energy method relies on the aerodynamical damping ratios (or equivalently, the aerodynamic work done per cycle) obtained in this way. A negative sign or the damping (or the aerodynamic work) indicates that the corresponding mode is prone to flutter vibrations. 11The transfer matrix contains also non-zero off-diagonal terms that describe the conservative and dissipative cross-coupling among the modes. This cross-coupling does not have to be reciprocal such that the transfer matrix is generally not symmetric.

Numerical computation As in the case of external aeroelastic forces, the aeroelastic coupling is determined using computational fluid dynamics. To this end, a set of modes is first defined for each IBPA to be considered in the coupling formulation. Here, the term 'modes' has a rather general meaning and is not necessarily restricted to natural modes of vibration. Instead, a set of representative modes can be selected, which comprises general vibrational deflection shapes that are deemed relevant for the expected vibration behavior. These modes do not have to coincide with the component modes used for the reduced description of the structural dynamics. However, the representative modes should ideally live in the subspace spanned by the component modes assembled in T , such that the aeroelastic coupling can be well-described in this basis.

In the second step, an unsteady fluid dynamics computation is carried out, see e. g. [105,[START_REF] Kersken | Time-Linearized and Time-Accurate 3D RANS Methods for Aeroelastic Analysis in Turbomachinery[END_REF]. The dynamical displacement of the blades is prescribed, separately for every considered mode, with oscillation frequency Ω.

The resulting fluid pressure is then integrated over the blade surface to determine the consistent nodal forces similar to Eq. [START_REF] Charleux | Threedimensional Multiharmonic Analysis of Contact and Friction in Dovetail Joints[END_REF]. However, in this case, it is important to account for the time variation of the surface normal, so that η S (x, t) is no longer constant but depends on time t. Moreover, only the fundamental component of the dynamic forces is usually considered, which corresponds to the imposed traveling wave motion. Finally, the forces are projected onto the set of representative modes to determine the aeroelastic transfer matrix tw G ai,k for each relevant IBPA θ k . If applicable, an according coordinate transform between representative and component modes needs to be taken into account.

The transfer matrix tw G ai,k depends on the oscillation frequency Ω used in the above described analysis. This frequency-dependence is generally nonlinear in nature, but it is sometimes linearized if the considered frequency range is comparatively small [START_REF] He | Fluid-structural coupling effects on the dynamics of mistuned bladed disks[END_REF]. It should be emphasized that the oscillation frequency is generally independent of the rotational speed Ω rot , and does not have to coincide with a particular natural frequency.

In the framework of the high-order harmonic balance method, presented in the next section, higher order frequencies nΩ are taken into account besides the fundamental oscillation frequency Ω. This should be considered in the aeroelastic coupling definition.

Comment on integrated approaches

It should be noted that integrated approaches exist for the closely coupled analysis of fluid-structure interactions [101, [START_REF] Dowell | Modeling of fluid-structure interaction[END_REF][START_REF] Kamakoti | Fluid-structure interaction for aeroelastic applications[END_REF]. However, these approaches have currently a limited applicability in the field of nonlinearlycoupled bladed disks, due to their prohibitive computational effort. Hence, these approaches are not presented in this work. The strategies for computing aerodynamical external forces and linearized aeroelastic coupling forces are presented in the following two subsubsections, respectively.

Simulation

In the previous section, a model for the structural dynamic behavior of bladed disks was presented, which considers the nonlinear contact interactions in the friction joints as well as the aeroelastic effects introduced by the surrounding fluid flow. In the present section, simulation approaches are presented. The purpose of these approaches is to determine the vibration behavior, i. e. , the time evolution of the coordinates based on presented model. The approaches can be applied to the equations of motion in physical coordinates, Eq. ( 2), in traveling wave coordinates, Eq. ( 4), and their projected variants (limited to the subspace of the component modes). To keep the formulations general, the vector u is used in this section, which shall represent the respective set of coordinates. The specific features in the case of traveling wave coordinates are highlighted, where applicable. This section is organized as follows.

In Subsection 3.1, the dynamic regime of interest is defined, for which the vibration behavior is to be determined. This definition is of crucial importance for the choice of suitable simulation techniques. In Subsection 3.2, an overview is given on general simulation methods capable of computing both non-periodic as well as periodic vibration regimes. Fourier methods represent a particularly popular family of methods capable of directly computing steady-state vibrations; they are addressed in detail in Subsection 3.3. These methods give rise to nonlinear equations. Appropriate numerical methods for the computation of solutions, their continuation and the calculation of their branching behavior are presented in Subsection 3.4. Finally, advanced simulation techniques of particular relevance to frictiondamped bladed disks are addressed in Subsection 3.5.

Dynamic regimes of interest

In this subsection, we describe the vibration regimes of particular interest in the design of bladed disks. To this end, vibration phenomena are categorized in Subsubsection 3.1.1. The category of periodic, traveling wavetype vibrations is of foremost importance, and it is described in detail in Subsubsection 3.1.2.

Vibration mechanisms and their dynamical features

In general, vibrations can be categorized according to the nature of the excitation and the character of the vibration response, see Fig. 10a and b respectively. Vibrations of bladed disks can be induced by either sustained or event-like excitation mechanisms. An example for an event-like excitation mechanism is the impulsive impingement of foreign or domestic objects on the blades. Also, the rubbing between blades and the casing is normally considered to be of event-like nature.

Event-like excitations typically lead to free vibrations and are not further discussed in this work. On the other hand, the important aeroelastic excitation mechanisms described in Subsection 2.4 represent sustained excitation mechanisms. Sustained excitation mechanisms can lead to forced vibrations, self-excited vibrations or parameter-excited vibrations. The most common source of forced vibrations is the engine order excitation introduced in Subsection 2.4. Here, the vibration behavior in the neighborhood of a specific resonance is typically of interest, i. e. , near a particular coincidence of an excitation frequency mΩ rot with one of the system's natural frequencies. The most common source of self-excited vibrations are aeroelastic instabilities which lead to flutter vibrations. Parametric excitation can be caused by the rotation in the gravity field, and certain other rotordynamical phenomena. Parameter-excited vibrations currently play a minor role in the design of bladed disks.

In the case of forced vibrations, the excitation can be further grouped into stationary and transient excitation. An example of an transient excitation is the external forcing during run-up or run-down of the machine. While, in this case, the excitation is typically regarded as sustained, its frequency spectrum varies with the rotational speed. From the design perspective, stationary excitation forces are particularly critical, since the vibrations can fully develop to a possibly high level and persist, which can cause fatigue. In addition, situations of resonance passages are of interest, where the timevarying excitation frequency coincides with one of the system's natural frequencies. While the vibrations are a priori transient, excessive vibration levels could still be reached, if the resonance is passed sufficiently slowly.

Resonance passages are attested an increasingly important role in the design of industrial turbines for power generation, since more run-ups and run-downs are required in the course of the transition to more renewable energies [START_REF] Bonhage | Transient amplitude amplification of mistuned blisks[END_REF]. The by far most important form of external forces are periodic excitation forces. We will focus on periodic excitation forces in this work, as these are of primal relevance in the design of bladed disks with regard to high cycle fatigue, and these represent one of the main motivations for the use of friction damping.

Certain aerodynamical phenomena may lead to random excitation, i. e. , stationary, non-periodic excitation where the external forces are described in terms of their stochastic characteristics rather than their exact timedependence (e. g. Gaussian noise).

In the presence of periodic external forces or aeroelastic instabilities, one can further distinguish between transient and steady-state vibrations, cf. Fig. 10b. Transient vibrations refer to the transitional regime from certain initial conditions until a steady state is reached or the vibrations grow unboundedly. Steady-state vibrations are typically of most importance. They can be either periodic or non-periodic, namely quasi-periodic or chaotic. Periodic vibrations are certainly the most important case. In bladed disks, they occur in particular as forced vibrations or flutter-induced limit cycle oscillations. In the case of flutter, it should be remarked that it is not uncommon to have multiple unstable traveling wave forms within a specific mode family. The resulting nonlinear vibration behavior in this case is still a lively area of research. First investigations suggest that multiple distinct steady states are possible in this case [102,103]. Under the condition of a rotationally periodic structure, each of these steady states corresponds to a limit cycle oscillation, where the vibration form is dominated by one of the unstable wave numbers. The basins of attraction of the limit cycles associated with the most unstable wave numbers tend to be the largest ones.

In the special situation where flutter occurs in the neighborhood of a possible external resonance, both excitation sources should be considered. First investigations suggest that there exists a lock-in frequency region close to the resonance coincidence point [START_REF] Corral | Conceptual Analysis of the Non-Linear Forced Response of Aerodynamically Unstable Bladed-Discs[END_REF]. In this lock-in region, the vibration is synchronous to the external forcing, and thus oscillates with an integer multiple of the rotational speed. In a certain distance from resonance, the autonomous behavior prevails, and the vibration is no longer synchronous to the rotor speed. According to [START_REF] Corral | Conceptual Analysis of the Non-Linear Forced Response of Aerodynamically Unstable Bladed-Discs[END_REF], the vibration behavior is periodic in a wide frequency range, with a potential chaotic transition regime In summary, the by far most important dynamic regime of bladed disks is that of periodic vibrations. Periodic vibrations are relevant both in the case of forced response and flutter. In the next subsubsection, we will discuss a particularly important form of periodic vibrations of bladed disks, namely traveling wave-type vibrations.

Periodic traveling wave-type vibrations

Consider a periodic vibration with fundamental frequency Ω,

u(t) = { ∞ n=0 U n e inΩt } . (19)
where U n are complex-valued amplitudes. A particularly important special case of periodic vibrations in rotationally periodic structures are traveling wave-type vibrations. The general notion of traveling waves was already introduced in A. A periodic vibration is considered to be of traveling wave-type, if the generalized coordinates (l) u of sector l exhibit a constant time lag l∆t with respect to the coordinates (0) u of the reference sector,

(l) u(t) = (0) u(t + l∆t) l ∈ [0, n s -1] . (20)
Herein, ∆t is the constant inter-sector time lag. Eq. ( 20) shall be valid for any choice of the reference sector. This implies that (0) u(t + n s ∆t) = (0) u(t). Therefore, ∆t can only assume discrete values ∆t = m 2π Ωns with an arbitrary integer m ∈ Z. Due to aliasing, it is sufficient to consider only wave forms corresponding to integers within the interval [0, n s -1]. We denote the integer Fig. 11: Illustration of the phase lag boundary condition defined in Eq. ( 21) m 0 ∈ [0, n s -1] fundamental wave number. By substituting Eq. ( 20) into Eq. ( 19), a relationship between the complex amplitude vector (l) U n of sector l and that of the reference sector can be established,

(l) U n = (0) U n e i 2π ns m0nl . ( 21
)
This relationship is illustrated in Fig. 11. This way, Eq. ( 20) can be understood as a constraint that only permits specific wave forms, namely those of traveling type. To provide more insight into this constraint, consider the following extreme cases:

1. A single blade vibrates, while all other blades remain static for all times. This vibration behavior is not captured by Eq. ( 20). Indeed, it represents a localized standing wave (soliton) from the perspective of the bladed disk.

2. A single blade oscillates, one after the other. This vibration behavior is in fact captured by Eq. ( 20). It represents a localized wave, traveling around the bladed disk. 3. All blades oscillate synchronously with the same magnitude. This is a standing wave in the rotating frame of reference. However it satisfies Eq. ( 20) since it can also be represented as the special case of a traveling wave with m 0 = 0 (due to the aliasing effect).

Traveling wave-type vibrations are expected if (a) the mechanical structure is rotationally periodic, (b) the aeroelastic effects do not introduce coupling among different traveling wave forms, and (c) the external forcing, if any, is of traveling wave type. The rotational periodicity of the mechanical structure, condition (a), implies that the effect of mistuning remains negligible 12 . Condition (b) implies that the aeroelastic coupling influences the traveling wave vibration forms individually. This is a common assumption, in particular, if the coupling is assumed to be of linear type, cf. Subsection 2.4. Finally, condition (c) is often satisfied since the most common form of external forcing is the engine order excitation, which takes the form of a traveling wave.

It must be noted that the above stated conditions are not sufficient, since strongly nonlinear effects can generally impede the emergence of traveling waves. As discussed in Subsubsection 2.1.3, however, the effects caused by the local contact interactions in friction-damped bladed disks are commonly assumed to remain weak or moderate. Hence, it is typically expected that under the above stated conditions, the traveling wave nature of system and forcing propagates to the vibration response. It is, however, taken into account that the local nonlinear effects may introduce certain higher-order components in the vibration behavior (in time and space).

The relation [START_REF] Claeys | Experiments and numerical simulations of nonlinear vibration responses of an assembly with friction joints -Application on a test structure named 'Harmony[END_REF] can be exploited to greatly simplify the simulation in the case of traveling wave-type vibrations. In accordance with the above stated assumption of rotational periodicity, the form of the contact laws and associated parameters are identical for each sector. Eq. ( 21) defines a phase shift between the input displacements for the nonlinear forces of each sector. This implies that the harmonics F c,n of the nonlinear forces strictly satisfy the same relation as the coordinates U n in Eq. ( 21). Hence, the nonlinear forces also take the form of periodic traveling waves. The nonlinear forces f c (t) acting on the global structure can thus be expressed in terms of the harmonic components (0) F c,n

of the nonlinear force acting on the reference sector,

f c (t) = ∞ n=0 √ n s w kn+1 ⊗ (0) F c,n e inΩt , (22)
where k n is defined by the congruence rule 13 ,

k n = nm 0 mod n s . (23)
Lastly, the harmonic components of the nonlinear force are transformed to traveling wave coordinates,

tw f c,j (t) = w H j+1 c(t) = n∈Kj √ n s (0) F c,n e inΩt . ( 24
)
with

K j = { m | mm 0 mod n s = j, m ∈ [0, ∞[ }.
Note that Eq. ( 24) only involves the contact forces at the interfaces of the reference sector, whereas in the general case of periodic or non-periodic vibrations, all contact interfaces of the bladed disk have to be considered individually. The problem of analyzing the vibrations of the global structure is thus reduced to the analysis of a reference sector. Therefore, the assumption of traveling wave-type vibrations is a very substantial one.

General methods for non-periodic and periodic motions

Exact solutions of the nonlinear equations of motion are limited to rather simplified special cases, which are of rather limited usefulness in the design of bladed disks coupled by mechanical joints 14 . Hence, only approximate solutions will be discussed in this work. In general, asymptotic solutions can be obtained analytically or semi-analytically by means of perturbation techniques. However, the application of such approaches to problems involving three-dimensional, possibly non-smooth contact interactions is a rather difficult issue. More importantly, these approaches have a comparatively limited range of applicability. Instead, numerical methods are usually employed for the vibration prediction of nonlinearly-coupled bladed disks.

The most versatile class of numerical simulation methods is the class of time integration methods. Starting 13 For enlightening illustrations of this rule, the reader is referred to [113] 14 It appears to be a common belief that the equations of motion can be solved exactly if piecewise linear contact laws are considered. Indeed, the set of ordinary differential equations becomes piecewise integrable. However, the transition times between the different contact states (stick, slip, liftoff) are generally not a priori known and need to be determined from the transition conditions. The latter are usually transcendental equations in the unknown transition times, rendering an exact solution impossible.

from a given set of initial values, the unknown time evolution of the generalized coordinates in accordance with the differential equations of motion is determined, successively, via appropriate quadrature rules until a specified stop time is reached. For contact problems, time integration methods can be grouped into event-driven and time stepping algorithms. Event-driven algorithms aim at directly resolving the transitions between different contact states. This becomes a comparatively expensive strategy if numerous events occur [START_REF] Feeny | Chaos in a Forced Dry-Friction Oscillator: Experiments and Numerical Modelling[END_REF]129]. On the other hand, event-driven algorithms lead to more relatively accurate results, since they do not suffer from the inherent discretization errors associated with time stepping algorithms. Time integration methods can also be categorized according to whether their quadrature rule can be given in explicit or only in implicit form.

Explicit methods are more prone to numerical stability problems, while implicit methods are known for inaccuracies associated with numerical damping [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF].

As discussed in Subsection 2.2, contact problems often lead to numerically stiff ODEs. Even if the contact problem is solved in its non-smooth formulation using appropriate methods, comparatively small time steps are required to ensure sufficient accuracy, which leads to considerable computational effort [131]. As an example, Phadke and Berger [130] simulated the steady-state forced response of a blade model with underplatform dampers by means of time step integration technique using a conventional finite element tool. To evaluate the vibration response for a single frequency point, a computation time of 40 hours was necessary (in 2008). This clearly demonstrates the prohibitive computational cost associated with time integration methods. Therefore, these methods cannot be applied in the context of comprehensive parameter studies, required for the design of bladed disks coupled by friction joints. However, time integration methods are widely used as reference for other approximate methods.

As discussed in Subsection 3.1, periodic vibrations are of primary importance in the design of bladed disks with mechanical joints. For their efficient computation, specific methods have been developed. The central idea of these approaches is the reformulation of the initial value problem into a two-point boundary value problem, where periodicity is enforced by requiring equality of the generalized coordinates and velocities at the beginning and at the end of the period. The methods for the solution of this boundary value problem can be grouped as follows:

(a) Shooting methods: The purpose of shooting methods is to find the appropriate initial values of the generalized coordinates and velocities which satisfy periodicity boundary condition. The problem is for-mulated as a set of nonlinear equations in the unknown initial values. For given initial values, the values at the end of the period are still determined by means of time integration. (b) Spectral methods: Here, the generalized coordinates and velocities are approximated in terms of suitably chosen base functions and associated coefficients. Periodic base functions are used, in order to satisfy the periodicity boundary condition a priori. Since the number of considered base functions is finite in practice, the equations of motion can in general not be satisfied at all times. Instead, the time-dependent error term is made orthogonal with respect to suitably selected weight (or test) functions. This projection gives rise to a set of nonlinear algebraic equations in the unknown coefficients.

In contrast to direct time integration from given initial values, the transient regime is not explicitly computed in the case of shooting methods. This is an important advantage, considering that the transient regime is comparatively long due to the typically weak damping [START_REF] Ferri | Frequency Domain Solutions to Multi-Degree-of-Freedom, Dry Friction Damped Systems[END_REF]. Owing to the use of time integration, shooting methods are applicable to various problems, including those involving non-smooth contact laws. It should be remarked, however, that shooting methods are only rarely applied to nonlinearly-coupled bladed disks [START_REF] Feeny | Chaos in a Forced Dry-Friction Oscillator: Experiments and Numerical Modelling[END_REF]. Many variants of spectral methods are available [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF].

Compared to shooting methods, they do in general not rely on time integration. If the base and weight functions are carefully selected, a high accuracy and a comparatively high efficiency can be achieved with spectral methods. An important category of spectral methods are Galerkin methods, where the base functions are used as weights. Another category are collocation methods, where the error term is enforced to vanish at certain time instants [START_REF] Hesthaven | Spectral methods for time-dependent problems[END_REF]. A particularly popular Galerkin method is the harmonic balance method, where harmonic base functions are used. This method is widely used in the vibration prediction of friction-damped bladed disks, and is presented in the next subsection.

The harmonic balance method

The harmonic balance method 15 is well suited for the computation of periodic solutions of ODEs [162, [START_REF] Cesari | Functional Analysis and Periodic Solutions of Nonlinear Differential Equations[END_REF]111].

As discussed in the previous subsection, the harmonic balance method is a Galerkin method with harmonic base functions. Hence, the generalized coordinates u(t) are expanded in a truncated Fourier series,

u(t) ≈ { n∈H U n e inΩt } . (25)
Herein, n are the (temporal) harmonic indices, Ω is the fundamental vibration frequency and U n are the complex-valued amplitudes (often simply referred to as harmonics). H is the set of considered harmonic components. In general, the size of the set H is referred to as the number of harmonics. In the simplest case, all harmonics up to a specified harmonic order n h are considered; i. e. , H = {0, . . . , n h }.

As mentioned before, the substitution of Eq. ( 25) into the equations of motion produces an error term. In accordance with the Galerkin idea, this error term should not have a component in the subspace spanned by the base functions. This is achieved by making the error term orthogonal to the base functions, i. e. , by ensuring that the Fourier components of the residual term are zero (Fourier-Galerkin projection). To this end, the different terms in the equations of motion are represented by truncated Fourier series. This gives rise to a set of nonlinear algebraic equations in the unknown Fourier coefficients U n ,

-(nΩ) 2 M n + inΩD n + K n + G ai,n (nΩ) Sn(Ω) U n + F c,n (U 0 , . . . , U n h , Ω) = F ae,n ∀ n ∈ H , (26)
Herein, S n is the so-called dynamic stiffness (or impedance) matrix with regard to the n-th harmonic. Regarding the structural matrices M n , D n , K n and the aeroelastic transfer matrix G ai,n one needs to distinguish between physical and traveling wave coordinates. In the case of physical coordinates, these matrices do not depend explicitly on n. In the case of traveling wave coordinates, the situation is slightly more complicated. The corresponding traveling wave structural matrices tw K k , tw D k and tw M k , as defined in B, have to be inserted. Thereby, the (spatial) wave number k and the temporal harmonic index n are related via the congruence rule k = nm 0 mod n s , cf. Eq. ( 23). Here, the fundamental wave number m 0 is defined as discussed in Subsubsection 3.1.2. Moreover, G ai,n in Eq. ( 26) is replaced by tw G ai,n , as defined in Eq. ( 18). In general, the complex amplitudes F c,n , F ae,n of the forces are related to their respective time-domain counterpart by the Fourier transform. However, F ae,n is readily available in the frequency domain, see Eq. ( 17). The treatment of the contact forces requires special attention and is addressed in Subsubsection 3.3.5.

The fundamental wave number m 0 and the fundamental oscillation frequency Ω depend on the considered load case. In the case of forced response to engine order excitation, typically the fundamental wave number is equal to the fundamental engine order and Ω = m 0 Ω rot . If period doubling or general period multiplication occur, m 0 and Ω have to be adjusted accordingly. In the case of flutter-induced limit-cycle vibrations, m 0 is equal to the considered (unstable) wave number, or, equivalently, the corresponding IBPA, cf. Tab. 3. The fundamental frequency of the limit cycle oscillation is not directly linked to the rotation speed and has to be treated as an unknown. As a consequence, the problem in Eqs. ( 26) is under-determined. In this autonomous case, the absolute phase of the oscillation is arbitrary. For normalization, a phase constraint can thus be imposed such that the number of equations matches the number of unknowns again. The Eqs. [START_REF] Corral | Conceptual Analysis of the Non-Linear Forced Response of Aerodynamically Unstable Bladed-Discs[END_REF] for each harmonic can be assembled into a global set of equations of the form,

S (Ω) U + F c (U , Ω) = F ae , (27)
or

U + H (Ω) F c (U , Ω) = U ae . (28)
Herein, S denotes the (global) dynamic stiffness matrix, and H denotes the dynamic compliance (or receptance, or FRF) matrix, which is related to S by HS = I. Eq. ( 28) is obtained from Eq. ( 27) by leftmultiplication with H. It should be noted that Eq. ( 27) and Eq. ( 28) are mathematically fully equivalent, as long as S is invertible. While Eq. ( 27) represents a balance of generalized forces, Eq. ( 28) balances generalized displacements, which can be interpreted as the requirement that the actual displacement should be compatible with the displacement induced by the external forces (U ae = HF ae ) as well as the displacement induced by the nonlinear forces (HF c ). The formulation in Eq. ( 28) is sometimes preferred if the FRF matrix H is expressed directly, e. g. in terms of the natural modes and possibly the static residual compliance [122]. Methods for the numerical solution of nonlinear equations ( 27) or [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF], including methods for the continuation of solution branches and computation of the branching behavior are presented in Subsection 3.4. Before that, we discuss different theoretical and computational aspects related to the harmonic balance method in the following three subsubsections.

Critique of the harmonic balance method

The approximation of u(t), defined in Eq. ( 25), is perfectly smooth, i. e. , the time-derivatives exist and are continuous for any order (differentiability class C ∞ ). The contact forces, on the other hand, have typically a limited degree of smoothness. Suppose that the contact forces are C 0 functions. As a consequence, the exact solution u(t) will be a C 2 function (ü(t) and u(t) are C 0 and C 1 , respectively). Hence, the exact solution does not live in the space spanned by the harmonic base functions, even if an infinite number of base functions was considered. Besides this mathematical contradiction, the harmonic balance method may suffer from poor convergence behavior. Specifically, it can be shown that the harmonic balance method converges only point-wise and at a comparatively low rate in the presence of non-smooth forces [START_REF] Kim | Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity[END_REF]. In the neighborhood of discontinuities, the approximation exhibits oscillatory behavior, known as the Gibbs phenomenon.

In spite of the above mentioned difficulties, the harmonic balance method is commonly applied to problems involving non-smooth contact forces. It should be noted that one is mainly interested in predicting the global vibration behavior rather than in the detailed resolution of the locally non-smooth contact interactions. Owing to the 'filtering' behavior of the elastic structure, the global deformation behavior exhibits a higher degree of smoothness than the local relative deformations, velocities and forces in the contact interface. Numerical indicate that the global vibration behavior can often be predicted with sufficient accuracy, as compared with time integration of the nonsmooth model [110]. This might explain the great success of the harmonic balance method in the field of jointed structures.

Choice of the harmonic order and selection of harmonics

The number of considered harmonics should be as large as necessary to approximate the exact solution with sufficient accuracy, cf. Fig. 12. In order to avoid spurious computational effort, on the other hand, the number of harmonics should be as small as possible. A minimum set of harmonics is dictated by the considered excitation. Those harmonics that are directly driven, or where a significant response is expected, should always be considered.

Example A non-trivial, but, at the same time, not uncommon example should be mentioned here: the external forcing with different harmonics. Suppose that the external forcing of engine order type is in accordance with Eq. ( 16), and only the engine orders m 1 = 20 and m 2 = 21 provide a significant loading; i. e. , M = Fig. 12: Example of the convergence of the harmonic balance method with respect to increasing harmonic order n h , with H = {0, . . . , n h } {m 1 , m 2 }. In this case, the fundamental engine order is m 0 = 1 (defined as the greatest common divisor).

It is certainly necessary to consider the 20th and the 21st harmonic, in order to capture the excitation. Besides multiples of 20 and 21, combination harmonics z 1 m 1 + z 2 m 2 with integers z 1 , z 2 ∈ Z become relevant in general. Determining the smallest set H that captures the dynamic behavior with satisfying accuracy is not trivial in such a situation.

Several strategies have been proposed to automatically adjust the harmonic order n h during the solution process, or, more generally, to adapt the set H of considered harmonics. A simple approach is to determine the approximation for an initial harmonic order and for a higher one. If the difference does not exceed a certain tolerance with respect to an appropriate measure, the initial harmonic order is sufficient, otherwise not. For instance, a measure based on the strain energy was proposed in [START_REF] Jaumouillé | An Adaptive Harmonic Balance Method for Predicting the Nonlinear Dynamic Responses of Mechanical Systems -Application to Bolted Structures[END_REF]. Grolet and Thouverez [START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF] suggested to estimate how the spectral energy distribution varies in the progress of the continuation of the solution. Specific harmonics can then be included or excluded depending on whether its associated energy is predicted to exceed or fall below a certain threshold, respectively. It can be useful to select different harmonic orders for different coordinates. If the nonlinearities are localized in the contact interfaces, it is usually necessary to approximate the associated coordinates with a comparatively high number of harmonics, while a much smaller number is often sufficient for the remaining coordinates.

Experiments and simulations of friction-damped systems often suggest that the fundamental harmonic component of the (global) response is predominant within the dynamic regime of interest [START_REF] Ji | Experimental Investigation of the Dynamic Characteristics of the Damped Blade[END_REF]. Hence, the singleterm variant, where only the harmonic n = 1 is retained in the balance [START_REF] Corral | Conceptual Analysis of the Non-Linear Forced Response of Aerodynamically Unstable Bladed-Discs[END_REF], has been very popular for a long time [140, [START_REF] Chen | Periodic Forced Response of Structures Having Three-Dimensional Frictional Constraints[END_REF]115,[START_REF] Borrajo | Analytical Formulation of the Jacobian Matrix for Non-Linear Calculation of the Forced Response of Turbine Blade Assemblies with Wedge Friction dampers[END_REF]. However, accounting for higher harmonics can considerably influence the fundamental harmonic component of the response [START_REF] Coudeyras | A New Treatment for Predicting the Self-Excited Vibrations of Nonlinear Systems with Frictional Interfaces: The Constrained Harmonic Balance Method, with application to disc brake squeal[END_REF]. The higher-order contributions are important for the accurate resolution of the local stick, slip, and lift-off behavior in the contact interface [START_REF] Chen | Periodic Forced Response of Structures Having Three-Dimensional Frictional Constraints[END_REF]. Furthermore, it is crucial to consider the static balance, i. e. , the zeroth harmonic in Eq. ( 26) for contact problems. Otherwise, the effects of settling and realignment of the contact interface [START_REF] Firrone | Underplatform Dampers for Turbine Blades: The Effect of Damper Static Balance on the Blade Dynamics[END_REF][START_REF] Firrone | The Effect of Underplatform Dampers on the Forced Response of Bladed Disks by a Coupled Static/Dynamic Harmonic Balance Method[END_REF]180,139] and the static deflection typical for breathing contact which involves dynamic opening and closing of the contact. Hence, the single-term variant can suffer from poor accuracy [START_REF] Ferri | Frequency Domain Solutions to Multi-Degree-of-Freedom, Dry Friction Damped Systems[END_REF]175]. Particularly for the case of small clearances or low normal pressures [175], so-called subor superharmonic resonances are reported for numerical as well as for experimental investigations [170, [START_REF] Berthillier | Blades Forced Response Analysis with Friction Dampers[END_REF][START_REF] Chen | Periodic Forced Response of Structures Having Three-Dimensional Frictional Constraints[END_REF][START_REF] Chen | Periodic Response of Blades Having Three-Dimensional Nonlinear Shroud Constraints[END_REF][START_REF] Krack | Multiharmonic Analysis and Design of Shroud Friction Joints of Bladed Disks Subject to Microslip[END_REF].

In these dynamic regimes, several modes and multiple frequencies (not only the fundamental one) contribute substantially to the global vibration behavior. In summary, the appropriate selection of harmonics depends on the considered dynamic regime, the contact formulation (see discussion in Subsubsection 3.3.5), and the desired accuracy. This makes it impossible to provide a general recommendation for all cases.

Stability of periodic motions

The periodic motion computed by means of the harmonic balance method is not necessarily stable. In this work, stability refers to the conventional concept of asymptotic stability of a periodic orbit in the phase space [160]. If the periodic orbit is stable, all trajectories starting in a sufficiently small neighborhood of the orbit will converge towards it. Otherwise, in the presence of slightest perturbations, a motion initiated on this orbit will diverge, and eventually approach a stable limit state in another region in phase space. Hence, the stability indicates whether the vibration behavior can be expected in the real world, where small perturbations are inevitable. It is therefore relevant to assess the stability of the computed periodic motion. The Floquet theorem can be utilized to assess a posteriori whether a found periodic motion is asymptotically stable. Two different computational methods for stability analysis are highlighted at this point:

(a) Monodromy matrix based method : The stability can be inferred from the eigenvalues of the so-called monodromy matrix (Floquet multipliers). The monodromy matrix represents the mapping of infinitesimal perturbations at the beginning of the period to the end of the period. If there is an eigenvalue with a magnitude greater than unity, this indicates a perturbation would be magnified, and, hence, the motion is unstable. If, on the other hand, all eigenvalues are located within the unit disk in the complex plane, the motion is asymptotically stable. (b) Hill method : The departure point for the Hill method is the linearization of the dynamic behavior around the periodic motion. This leads to a linear ordinary differential equation with periodic coefficient matrices. Using the Hill method, the problem can be transformed to a quadratic eigenvalue problem in the frequency domain [START_REF] Groll | The Harmonic Balance Method With Arc-Length Continuation in Rotor/Stator Contact Problems[END_REF]. The eigenvalues are the so-called Floquet exponents. If all eigenvalues have negative real part, the periodic motion is stable, whereas it is unstable if any of the eigenvalues exhibits a positive real part.

For method (a), the central task is the computation of the monodromy matrix. This usually involves numerical time integrations. Hence, this approach is comparatively expensive. On the other hand, this method is particularly popular in combination with the shooting method, where the monodromy matrix is often readily available from the calculation of the Jacobian of the residual [114,157,158]. The Hill method is well-suited for the frequency domain framework [START_REF] Groll | The Harmonic Balance Method With Arc-Length Continuation in Rotor/Stator Contact Problems[END_REF][START_REF] Kim | Super-and subharmonic response calculations for a torsional system with clearance nonlinearity using the harmonicbalanc e method[END_REF][START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. The coefficient matrices of the quadratic eigenvalue problem only depend on quantities readily available from the solution of the harmonic balance equations. A difficulty with the Hill method is the size of the quadratic eigenvalue problem, which scales with the number of harmonics. For large problems, the stability analysis based on the Hill method can become numerically more expensive than the computation of the periodic motion itself. Moreover, a much larger number of harmonics is commonly required for the accurate stability analysis compared to the computation of the vibration behavior [118]. Finally, a numerical difficulty is associated with the finite harmonic truncation of the Hill matrix [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. The number of eigenvalues obtained by the Hill method exceeds the number of Floquet multipliers by a factor corresponding to the number of harmonics. Owing to finite truncation, however, these eigenvalues are not identical, and some of them are more accurate than others. It is suggested to filter the eigenvalue spectrum in order to obtain only the most converged values [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF].

Under the variation of a system parameter, the periodic motion may change its stability, which is referred to as bifurcation. In the case of an unstable periodic motion, one or more new stable motions may emerge beyond a bifurcation point, each being either periodic or non-periodic. The different possible steady states of a nonlinear dynamical system are illustrated in the so- 13b. A Poincaré map illustrates the returns of a trajectory to a section through the phase space. A common example for bifurcation points are the turning points in the frequency response of nonlinear mechanical systems. Consider the case of a structure with opening contact; this system exhibits a softening-type nonlinear behavior. In this case, the amplitude-frequency curves of periodic forced vibrations near a specific resonance are typically bent towards the left, as indicated in Fig. 13a. As a consequence, there is a range with (at least) three possible amplitudes for the same frequency. Not all of the solutions are stable. At the turning points, saddle node bifurcations take place. The periodic motions associated with the overhanging branch, connecting the turning points, is typically unstable. In general, bifurcations may also give rise to emanating solution branches, and detached solution branches may be present, as illustrated in Fig. 13a. It can generally be stated that the analysis of the stability and the bifurcation behavior of periodic solutions, as well as the consideration of quasi-periodic and nonperiodic motions, are only rarely addressed in the field of nonlinearly-coupled bladed disks. In Subsection 3.4, methods are presented for the continuation and the computation of the branching behavior of solutions.

Extension to quasi-periodic motions

Besides periodic motions, chaotic and quasi-periodic motions are possible forms of steady-state vibrations of bladed disks. Quasi-periodic motions represent motions with multiple rationally independent (or noncommensurable) base frequencies. From a topological point of view, the trajectory of a quasi-periodic motion covers an invariant torus in the phase space, as opposed to the closed orbit in the case of a periodic motion.

In the presence of multi-frequency excitation with incommensurable frequencies, the vibration response is a priori quasi-periodic. But even in the presence of harmonic external forcing, nonlinearity can lead to the loss of stability of periodic motions, and give rise to quasiperiodic motions. Besides the excitation frequency, a new frequency emanates in this case, which contributes to the vibration behavior.

To compute quasi-periodic motions, both the shooting method and the harmonic balance method can be extended, see [START_REF] Ling | Quasi-periodic solutions calculated with the simple shooting technique[END_REF] and [142], respectively. In the case of Fourier methods, this extension is known as the multidimensional harmonic balance method (or variablecoefficient harmonic balance method). The central idea is to consider the rationally independent frequencies as fundamental frequencies. Hence, Ω and n in Eq. ( 25) become vectors of according dimensions, and nΩ is to be understood as inner product. Thus, besides multiples of either fundamental frequency, combination frequencies can be taken into account. In the case of traveling wave-type vibrations, the fundamental wave numbers associated with the fundamental frequencies are collected in the vector m 0 , and the congruence rule in Eq. ( 23) governs the corresponding wave numbers associated with hyper-index n, cf. e. g. [START_REF] Laxalde | Non-linear vibrations of multi-stage bladed disks systems with friction ring dampers[END_REF]. Accordingly, the one-dimensional (inverse) Fourier transform has to be replaced by its multi-dimensional counterpart, thus, acting on a multi-dimensional hyper-time. While hysteresis can be handled in the one-dimensional case, as will be detailed in Subsubsection 3.3.5, it represents a principal problem in hyper-time. Therefore, hysteretic nonlinearities, such as dry friction, have to be re-formu-lated by incorporating appropriate internal variables into the set of differential equations, such that the phase space does no longer have a 'memory' [START_REF] Krack | On the interaction of multiple traveling wave modes in the flutter vibrations of friction-damped tuned bladed disks[END_REF]. Finally, it should be mentioned that for assessing the asymptotic stability of quasi-periodic motions, the Floquet stability theorem can be accordingly extended [START_REF] Guskov | Harmonic Balance-Based Approach for Quasi-Periodic Motions and Stability Analysis[END_REF].

General strategies for the treatment of the nonlinear forces

Most harmonic balance formulations in the field of bladed disks with contact interfaces do not treat the nonlinear forces explicitly as unknowns. Hence, the nonlinear equations are solved only for the harmonics of the generalized coordinates U , as indicated in Eqs. ( 27)- [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF]. To this end, the harmonics of the nonlinear forces F c need to be expressed as a function of U . This is analogous to the primal assembly approach in dynamic substructuring, where the a priori unknown coupling forces are expressed as a function of the displacements, and not treated explicitly as unknowns. In this sense, a dual assembly is also possible, which means that the nonlinear forces F c are considered as additional unknowns.

Alternatively, the problem in Eqs. ( 27)-( 28) could be re-formulated in such a way that the nonlinear forces F c are treated as the only unknowns. It should be kept in mind that the nonlinear forces are typically of lower degree of smoothness than the generalized coordinates. Therefore, a larger number of harmonics would is typically necessary to approximate the nonlinear forces with satisfying accuracy. By treating the nonlinear forces implicity, this problem is circumvented, and a comparatively small number of harmonics is usually sufficient in the simulation.

In most cases, it is not possible to express the relation between the harmonics of the generalized coordinates U and the harmonics of the nonlinear forces F c in closed form. Instead, numerical procedures have to be utilized to evaluate the term F c (U ). This is a crucial task and often represents the bottleneck of the simulation procedure. A major problem is the case of set-valued relations between u and f c . This problem can be tackled by means of the Dynamic Lagrangian method presented in Subsubsection 3.3.6. Two aspects further complicate the computation of F c (U ): (a) The need to account for the different contact states (stick, slip, liftoff), and (b) the hysteretic character of the friction force. Regarding (a), the approaches can be grouped into event-driven and time stepping schemes, analogously to time integration methods. These approaches are presented in the following two paragraphs. The specific treatment of the hysteretic character is addressed in the third paragraph.

For convenience, we assume a regular relation of the form f c [u] in the following; the presented approaches apply accordingly in conjunction with the Dynamic Lagrangian method.

Event-driven schemes In the case of event-driven schemes, the key idea is to directly determine the transitions between the different contact states [124,[START_REF] Borrajo | Analytical Formulation of the Jacobian Matrix for Non-Linear Calculation of the Forced Response of Turbine Blade Assemblies with Wedge Friction dampers[END_REF][START_REF] Krack | A High-Order Harmonic Balance Method for Systems With Distinct States[END_REF]. This step involves the root finding of (functions of) trigonometric polynomials. In general, the roots of trigonometric polynomials are governed by transcendental equations, so that the computation of the transitions has to be carried out numerically. Once the transition time instants are known, the harmonics of the nonlinear forces F c can be computed by piecewise integration. As opposed to the computation of the transition times, this integration step can be carried out analytically for many systems. For the class of piecewise polynomial systems, this has been shown in [START_REF] Krack | A High-Order Harmonic Balance Method for Systems With Distinct States[END_REF] 16 . Note that piecewise polynomial systems include systems with unilateral elastic contact and one-dimensional elastic Coulomb friction, which are piecewise linear. However, the comparatively common case of spatial Coulomb friction, where the friction force needs to remain within the friction cone, is not included. In fact, an analytical integration is impossible in the case of spatial friction. Hence, spatial friction can only be tackled by considering the two coordinates, that span the tangential contact plane, as independent of each other, which represents a major limitation 17 . The analytical formulation results in a comparatively accurate description of the functional value F c and its derivatives of different order. This can be an important aspect in the context of design studies and sensitivity analyses. However, in accordance with the discussion in Subsection 3.2, event-driven schemes become comparatively expensive if numerous events occur (per vibration period). This is why they are less frequently applied to structures with a fine discretization of the contact area.

As described above, the Gibbs phenomenon leads to artificial oscillations of the Fourier approximation near discontinuities such as contact state transitions. This, in turn, may lead to artificial contact state transitions in addition to the (physical) contact state transitions. Hence, the character of the approximation, by design, gives rise to a large number of events, which makes event-driven schemes less attractive in the context of the Harmonic Balance method.

Fig. 14: Alternating-frequency-time scheme

Alternating frequency-time scheme The evaluation of F c using time stepping is referred to as alternating frequency-time (AFT) scheme [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF]. It relies on the evaluation of the nonlinear forces at discrete time steps and the conversion between time and frequency domain, as illustrated in Fig. 14. The evaluation of the nonlinear force laws in the time domain is usually straightforward. Different events are handled simply by piecewise definition of the nonlinear forces. The AFT scheme is very popular also beyond the realm of contact nonlinearities. Mathematically, the AFT scheme can be expressed as

F c = DFT [ f c (iDFT [U 0 , . . . , U n h ]) ] , (29)
where (i)DFT denotes the (inverse) discrete Fourier transform. It is particularly efficient to utilize the (inverse) Fast Fourier Transform ((i)FFT) for this step [START_REF] Cardona | A Multiharmonic Method for Non-Linear Vibration Analysis[END_REF].

Remark Note that for standard (i)FFT subroutines, a complete set of harmonics is usually required. For nontrivial sets H, the remaining harmonics U n of the generalized coordinates are set to zero. Accordingly, from the harmonics F c of the nonlinear forces obtained by the FFT, only those associated with the set H are considered in the nonlinear equations ( 27)- [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF].

Compared to the event-driven scheme, the AFT scheme exhibits an inherent discretization error. The accuracy and computational effort of the AFT scheme largely depend on the number of time steps per period. To avoid aliasing, the theoretical lower limit of this number is given by 2n h +1. To accurately resolve the contact state transitions, usually a much larger sampling frequency is needed.

Treatment of hysteresis If the nonlinear forces describe hysteretic effects, as in the case of dry friction, the force f c cannot be expressed as an explicit function of the generalized coordinates u and velocities u. Instead, the force is usually governed by a differential law such as the elastic Coulomb law defined in Eq. ( 11). Since periodic, steady-state vibrations are sought, one is only interested in the steady-state hysteresis cycle. A common strategy to determine the steady-state hysteresis cycle for given periodic input u(t), is to start from a certain point on the initial loading curve and to let f c evolve until reaching the steady-state hysteresis. For the AFT scheme, this can be achieved by evaluating Eq. ( 29) for several periods, and checking for periodicity of the forces (and/or further internal variables); in the case of elastic Coulomb friction, usually two periods are sufficient. For the event-driven scheme, state transitions are successively determined until the sequence of transitions starts to repeat itself.

Formulation of common contact formulations in the frequency domain

As stated in Subsection 2.2, the most popular contact formulations for bladed disks with mechanical joints are based on the Coulomb law for the frictional behavior, and the unilateral law for the normal contact. Specifically, two different formulations are commonly used: The contact laws are illustrated in Fig. 7. For a general discussion on these two formulations, we refer to Subsection 2.2. In the present subsubsection, it is shown how these formulations can be applied in the frequency domain.

In the case of the elastic formulation (a), the relation f c [u] is regular, such that the expression F c (U ) in Eqs. ( 27)-( 28) can be ad hoc evaluated by means of the general methods presented in the previous subsubsection. In contrast, the problem becomes a constrained optimization problem due to the set-valued force laws in case (b). Hence, a solution U is sought that satisfies the dynamic force equilibrium given in Eqs. ( 27)-( 28), as well as the contact constraints. This generally calls for specific algorithms well-suited for such problems. An appropriate method is the Dynamic Lagrangian (DL) method [110]. This method is similar in spirit as the well-known augmented Lagrangian method, but is specifically adapted to Fourier methods. An important difference to the augmented Lagrangian method is that the DL method is a primal method in the sense that the contact forces can be expressed as a function only of the primary displacement variables.

In the following, we present a unified formulation of both the elastic and the (DL based) rigid contact formulations. The difference between these formulations manifests itself merely in the definition of the predicted sticking force, which is the only input variable of the force computation algorithm. In Tab. 2, this predicted sticking force is defined in terms of its harmonic components for the different contact formulations. The force is given separately for the normal contact, Λ st,n , and the tangential contact, Λ st,t . The indices n and t refer to the restriction to the normal direction and the tangential directions, respectively, of a specific contact point. The fully three-dimensional case is considered here, where the normal contact influences the frictional behavior, and spatial friction is taken into account. Note that special cases such as frictional contact with constant normal load, frictionless unilateral contact, or frictional contact with decoupled tangential directions can be derived in a straight-forward manner. Contact is formulated on the force level, as opposed to the pressure level. Forces and pressures at contact point k are related by the associated area ∆A k , cf. Eq. ( 7). In the case of the elastic formulation, two different cases should be distinguished: (a) the coupling by a continuous elastic elastic layer, and (b) the coupling by bulk springs with stiffness values k n and k t . In case (a), the area-specific stiffness of the elastic layer (measured in N/m 3 ) has to be weighted with ∆A k to obtain the stiffness values k n and k t used in this subsubsection. In case (b), k n and k t are directly specified.

In the case of the DL formulation, penalty coefficients DL,n and DL,t are introduced. These penalty coefficients have a crucial influence on the mathematical conditioning, and, thus, the convergence behavior, but they do not have an immediate effect on the converged results. This is the substantial difference between the penalty coefficients needed in the DL formulation, compared to the stiffness coefficients involved in the elastic formulation, which generally has a distinct influence. In [START_REF] Charleux | Numerical and Experimental Study of Friction Damping in Blade Attachments of Rotating Bladed Disks[END_REF], it is recommended to use a value in the order of magnitude of the spectral radius of the Schur complement of the dynamic stiffness matrix. The DL formulation utilizes the term R lin , which represents the linear part of the residual of Eq. ( 27),

R lin := SU -F ae . (30)
Therefore, the non-contact forces, i. e. , elastic, inertia, damping and excitation forces, directly enter the contact force computation in the case of the DL formulation. Consider the case of fully sticking contact. Accord-ing to the DL formulation, g n = 0 and g t = 0 in this case. From this, it follows that Λ st,n = -R lin,n and Λ st,t = -R lin,t , which could also have been obtained by a direct elimination of the linear constraints. The consideration of non-contact forces in the contact force computation is an essential difference to the elastic formulation, where only the contact deformation enters the contact force computation. As the dynamic force balance given in Eq. ( 27) must hold at the solution point, the converged contact forces, of course, implicitly depend on the remaining forces, as in the case of the elastic formulation.

For frictional contact, a velocity-based variant can be used in the case of the DL formulation, as alternative to the displacement-based one, see the third row of Tab. 2. This has no counterpart in the elastic formulation. The matrix ∇ in Tab. 2 takes care of the time differentiation in the frequency domain; i. e. , ∇g t corresponds to the harmonic components of the tangential velocity.

Table 2: Definition of the predicted sticking force for different contact formulations

sticking force elastic Dynamic Lagrangian Λ st,n k n g n -R lin,n + DL,n g n Λ st,t (displacement-based) k t g t -R lin,t + DL,t g t Λ st,t (velocity-based) - -R lin,t + DL,t ∇g t
We exemplify the unified contact treatment in the framework of the AFT scheme. We would like to emphasize, however, that analogous formulations are possible with event-driven schemes. To this end, the frequencydomain quantities Λ st,n , and Λ st,t , are transformed into the time domain via the iFFT. The resulting time samples at time instant i are denoted λ st,n,i , λ st,t,i . Based on these time samples of the predicted sticking forces, the actual contact forces are determined. This procedure is presented in the following two paragraphs for normal and tangential contact, respectively.

Normal contact For the normal contact, the (actual) contact force λ n,i reads,

λ n,i = 0 if N 0 + λ st,n,i < 0 (separation) λ st,n,i if N 0 + λ st,n,i ≥ 0 (contact) . (31)
Herein, N 0 is the normal preload 18 at the considered contact point. The case of an initial clearance g n,0 is also captured in this formalism. In this case, negative values N 0 have to be used in accordance with either N 0 = -k n g n,0 for the elastic formulation, or N 0 = -DL,n g n,0 for the DL formulation. It can be verified that a negative normal preload corresponds to an offset of the normal contact gap by g n,0 , cf. Tab. 2.

Tangential contact Once the normal force is known for a specific contact point, the tangential contact force can be computed. In the case of displacement-based formulations (elastic or DL formulation), the following applies,

λ t,i =          λpre,t,i λ t,i-1 + λ st,t,i -λ st,t,i-1 (I) µ |λ n,i | λpre,t,i λpre,t,i (II) 0 (III) (I) if λ pre,t,i < µ |N 0 + λ n,i | (sticking) (II) if λ pre,t,i ≥ µ |N 0 + λ n,i | (sliding) (III) if N 0 + λ n,i = 0 (separation) . (32)
In contrast to the normal contact force, the tangential contact force λ t,i at time instant i depends on its value at the previous time instant, λ t,i-1 , resulting in an iterative scheme. This reflects the hysteretic character of the model. Eq. ( 32) has to be evaluated iteratively until the steady-state hysteresis cycle is reached.

As alternative to the displacement-based formulation, a velocity-based formulation can be used in the case of the DL formulation,

λ t,i =      λ st,t,i (I) µ |N 0 + λ n,i | λst,t,i λst,t,i (II) 0 (III) (I) if λ st,t,i < µ |N 0 + λ n,i | (sticking contact) (II) if λ st,t,i ≥ µ |N 0 + λ n,i | (slipping contact) (III) if N 0 + λ n,i = 0 (liftoff) . (33)
The obvious advantage is that no iterative evaluation of the force law is required in this case. According to [START_REF] Laxalde | Vibration Control for Integrally Bladed Disks Using Friction Ring Dampers[END_REF], however, this approach is not useful if the mean relative displacement of the contact interface is not known in advance.

Once the (steady-state) time histories of the forces λ n,i , λ t,i are known, the DFT is applied to determine the required harmonic components of the contact forces. Finally, the global nonlinear force vector F c can be assembled considering all contact interfaces.

Remark In the case of the elastic formulation, it is common to formulate the contact laws in terms of the gaps rather than using the predicted sticking force. Also, in the case of the DL formulation, the contact force is usually split into a frequency-domain term accounting for the dynamic force equilibrium and a time-domain term penalizing the violation of the contact constraints. Only the latter term is typically evaluated in the time domain. It can be verified that Eqs. ( 31)-( 33) are indeed fully equivalent to the common elastic and DL formulations. The alternative formulation was presented here to indicate the opportunity for the unified implementation of both contact formulations (using the same subroutine). Also, the notion of the predicted sticking force is deemed helpful to gain further insight into the interpretation of the Dynamic Lagrangian formulation and its relation to the elastic formulation. Finally, it should be noted that Eqs. ( 31) and ( 33) represent proximal point operations with respect to the admissible set of contact forces, applied to the predicted sticking force. This highlights the resemblance with formulations used in time integration schemes for the algebraic inclusion problem associated with the non-smooth contact laws [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF].

Numerical solution of the nonlinear equations

The harmonic balance method gives rise to a set of nonlinear algebraic equations ( 27)- [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF]. As described in Subsection 3.2, such equations also arise in the general case of spectral methods as well as the shooting method.

Hence, the problem of predicting periodic motions of nonlinear dynamical systems often involves the solution of nonlinear equations. In general, the set of equations can be written as

R(X) = 0 , (34)
where X is the vector of unknowns and R is the residual vector function, which measures the error if X is not the solution. The efficient solution of Eq. ( 34) in the light of vibration problems is the purpose of the present subsection.

Remark Many solution methods require real arithmetic, at least for the vector of unknowns X. To this end, the complex-valued unknowns U in Eqs. ( 27)-( 28) can be split into real and imaginary parts, and gathered in the vector X. The same can be done for the residual, in order to obtain a real-valued vector function R.

Exact solutions of the nonlinear equations are only possible in very special cases, and are not further discussed in this work. Numerical methods are commonly used to compute approximate solutions. In general, a multitude of solutions may exist for the same set of equations. The numerical methods can be grouped into global methods, which are suited to compute all solutions, and local methods, which are suited to compute a single solution in the neighborhood of the initial guess. Most methods used for the vibration prediction of nonlinear structures are local methods, and we therefore focus on local methods in the following. Global methods will be only briefly discussed in the context of isolated solution branches, in Subsubsection 3.4.5.

For local methods, an initial guess X 0 must be provided 19 . This initial guess is usually not a solution and therefore produces an error R(X 0 ) = 0. The method then successively (and often successfully) minimizes this error by adjusting the current guess. Certainly the most popular family of methods in the field of structural dynamics are variants of the Newton and the quasi-Newton method, presented in the following subsubsection. An interesting alternative are the so-called pseudotime solvers. The idea of these methods is to re-formulate the algebraic equations as an initial value problem, ∂X ∂t * + R(X) = 0, by adding a term ∂X ∂t * , where t * is a pseudo time variable. The steady solution of this initial value problem is characterized by ∂X ∂t * , and thus identical to the sought solution of the algebraic equations. The initial value problem is solved by means of numerical integration schemes, starting from the specified initial guess. Pseudo-time solvers can exhibit superior performance as compared with many other methods, in particular if a large number of unknowns is involved. These methods are popular in the field of computational fluid dynamics and fluid-structure interaction problems [106, 104], see [138] for an application to structural dynamics.

The Newton method and its related variants

A very popular local method for the solution of nonlinear algebraic equations is the Newton method 20 . An iteration of the Newton method involves the solution of a set of linear equations,

∂R ∂X T X=X k (X k+1 -X k ) = -R(X k ) . (35)
Herein, k denotes the current iteration number. Based on the current guess X k , the solution of the subproblem in Eq. ( 35) yields the next (and hopefully better) guess X k+1 . Starting from an initial guess X 0 , successive iterations are computed until the norm of the residual, R (X k ) < is smaller than a specified tolerance in a suitable norm. The Newton method utilizes the gradient ∂R ∂X T of the residual vector with respect to the vector of unknowns, the so-called Jacobian matrix. It can thus be characterized as a gradient-based method. Gradient-based methods are known for their great convergence behavior in the neighborhood of the solution. Such methods are therefore particularly well-suited if a good initial guess is available. Disadvantages of the classical Newton method are the fact that convergence is not always ensured, and the relatively large computational effort associated with the need to evaluate and to factorize the Jacobian in each iteration. In the following paragraphs, several improvements addressing these drawbacks are discussed. Several extensions and improvements of the classical Newton method have been proposed. Jacobian-free-Newton-Krylov methods combine a Newton-like procedure with Krylov subspace methods for solving the correction equations Eq. ( 35) and are well-adapted to large nonlinear problems [START_REF] Knoll | Jacobian-free Newton-Krylov methods: A survey of approaches and applications[END_REF]. Different improvements aim at ensuring global convergence, e. g. by equipping the method with a line-search or a trust-region algorithm. Such globally convergent variants of the Newton method are available in conventional computing toolboxes and libraries.

Quasi-Newton methods

The idea of the so-called quasi-Newton methods is to only approximate the Jacobian or its inverse, as opposed to computing these quantities accurately in each iteration. A simple strategy is to update the Jacobian not in every iteration. If only an approximation is used, inferior convergence behavior is expected. Thus, more iterations are typically required to achieve acceptable accuracy. However, since the average computation time per iteration is decreased, the total computation time of the method might be reduced. Another approach consists in approximating the inverse of the Jacobian, which is actually required to determine the next correction of the vector of unknowns in Eq. [START_REF] Feeny | Chaos in a Forced Dry-Friction Oscillator: Experiments and Numerical Modelling[END_REF]. To this end, the current approximate of the inverse Jacobian can be updated in each iteration, based on the quantities readily available. It is thus neither necessary to evaluate nor to factorize the Jacobian in each iteration. A common example is the Broyden-Fletcher-Goldfarb-Shanno method, see e. g. [133].

Analytical gradients For gradient-based methods, the solution process can be significantly accelerated by calculating the Jacobian ∂R ∂X T analytically, simultaneously to the evaluation of the residual R. This can greatly reduce the required computational cost compared to the common finite difference approximations of the Jacobian. The critical part is the derivation of the nonlinear forces F c . In the framework of the AFT scheme, this can be achieved by utilizing the linearity of the FFT and IFFT operators in Eq. ( 29),

∂F c ∂x = FFT ∂f c (u, u) ∂x + ∂f c (u, u) ∂u T • IFFT ∂U ∂x + ∂f (u, u) ∂ uT • IFFT ∂∇U ∂x . (36)
Herein, x is an arbitrary real-valued variable (not necessarily x ∈ X). The first term in Eq. ( 36) accounts for the case where f c depends explicitly on x, e. g. if x is a friction coefficient or a normal preload. The second and third term account for the case x ∈ {U } or x ∈ {U }; i. e. , when x is an element of the real or the imaginary part of the harmonic components U . Moreover, the third term also captures the case x = Ω, since ∇ depends on Ω (more specifically, is proportional to Ω). As discussed in Subsubsection 3.3.6, it is interesting to re-formulate the contact forces f c in terms of the harmonics of the predicted contact sticking forces, which, in turn, depend on X. In this case, the chain rule must be applied accordingly. For detailed examples, we refer to [START_REF] Berthillier | Blades Forced Response Analysis with Friction Dampers[END_REF]150].

Exploiting the sparsity of the nonlinear terms Frictiondamped systems are characterized by localized nonlinearities; i. e. , the nonlinear forces are confined to certain coupling interfaces at the joints. Provided that appropriate local coordinates are used for the description of the nonlinear coupling interfaces, as discussed in Subsection 2.3, the nonlinear forces f c act and depend on only a subset c g of the generalized coordinates u; i. e. , c g ⊂ u. As a consequence, the nonlinear terms Λ in Eqs. ( 27)-( 28) and the associated gradients are sparse. This sparsity can be exploited by an exact condensation procedure. The result is a reduced set of equations with a smaller number of unknowns. This can greatly improve the computational efficiency of the solution process. The condensation procedure involves computing (a large portion of) the dynamic compliance matrix, and appropriate care must be taken to avoid an expensive matrix inversion in every iteration. Details on this procedure are given in E.

Continuation of solution branches

During the analysis and the design of nonlinearly-coupled bladed disks, it is commonly relevant to determine the vibration behavior in certain ranges of a parameter ρ,

R(X, ρ) = 0 , ρ ∈ [ρ min , ρ max ] . (37)
Herein, ρ min and ρ max denote the lower and upper parameter limits, respectively, in between which the solution X is of interest. In the following, we limit the discussion to the case of scalar, real-valued parameters ρ ∈ R. For example, in the case of forced response analyses, a typical parameter is the oscillation frequency ρ = Ω. A set of connected solutions X(ρ) is denoted solution branch. Numerical (path) continuation methods facilitate the iterative computation of solution branches by accounting for already known solution points. A central difficulty in solving Eq. ( 37) is the possible existence of multiple solutions for the same value of ρ, as illustrated in the conceptual bifurcation diagram in Fig. 13a. The singularities at turning points can be overcome by means of continuation. In contrast, the detection of branching points and the switching to another branch require separate approaches. Finally, isolated branches may be present which cannot be directly found by means of local analysis. In the following, we discuss two popular continuation methods and elaborate on how to overcome the aforementioned difficulties encountered during continuation. A good review of suitable methods for continuation and bifurcation analysis, and related computational aspects is given in [148].

Fig. 15: Numerical continuation using a tangent predictor step and an arc length constraint for the corrector step

Predictor-corrector method The by far most common continuation technique is the predictor-corrector method. It is illustrated in Fig. 15. Starting from a given solution point, the next point is predicted by going a certain distance ∆s 'forward' into the estimated direction of the branch. For example, this direction can be the tangent or the secant at the solution branch. The step size taken in the predictor step can be adapted automatically, in order to avoid spurious computational effort resulting from too small steps, and ensuring fast convergence by choosing not tool large steps. The predicted point is generally not located on the solution branch. A correction step is therefore needed. To this end, the methods described in Subsubsection 3.4.1 can be used. In this step, the parameter ρ is generally considered as unknown. This makes Eq. ( 37) under-determined. An additional constraint equation is therefore introduced. A popular constraint is to ensure that the next solution point has the specified distance ∆s (arc length) from the current solution point 21 . Another popular constraint forces the solution to lie on the hyperplane through the predicted point, orthogonal to the tangent. Yet another constraint simply fixes one unknown to its predicted value (local parametrization). The different prediction and correction strategies can be combined. A very popular combination is the tangent predictor step with the arc length parametrization. Note that the sequential (or natural) continuation strategy, where the problem in Eq. ( 34) is simply solved for a sequence of parameter values, represents a special case of the predictor-corrector method. In this case, the previous solution point is used as prediction, and a local parametrization with always the same variable ρ is used. However, turning points with respect to the parameter ρ cannot be overcome with this strategy.

Asymptotic numerical method The tangent step in the predictor-corrector method represents an evaluation of the first-order Taylor series expansion around the current solution point. The error of this prediction can be generally reduced by increasing the order of the power series expansion. This is the central idea of the asymptotic numerical method (ANM) [START_REF] Cochelin | A High Order Purely Frequency-Based Harmonic Balance Formulation for Continuation of Periodic Solutions[END_REF]. In this method, typical expansion orders are in the range of 20 to 40. To formulate the expansions, derivatives of R with respect to Y T = X T ρ of according orders are required. To determine these derivatives, even for large orders, the ANM makes use of the special properties of quadraticorder multivariate polynomial equations. To this end, the equations need to be expressed in quadratic form

R(Y ) = R 0 + R 1 (Y ) + R 2 (Y , Y ) , (38)
where R 0 , R 1 and R 2 are constant, linear and quadratic functions in Y , respectively. Based on Eq. ( 38), the derivatives can be expressed analytically up to the desired order involving a recursive formula. The path can be expressed as a power series expansion in the arc length s. One advantage of the power series framework is that error estimations can be easily made. It is thus possible to determine the radius of convergence in terms 21 This method is sometimes referred to as the Riks method.

of a given tolerance. Therefore, the step size can be automatically adjusted in such a way that correction steps are no longer necessary. Moreover, by analyzing the power series, one can detect bifurcation points and determine the tangents to the emerging branches [START_REF] Cochelin | Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations[END_REF]. The present limitation to quadratic order problems [START_REF] Firrone | Underplatform Dampers for Turbine Blades: The Effect of Damper Static Balance on the Blade Dynamics[END_REF] is the main caveat of the ANM. It impedes the ad hoc application to systems involving dry friction and unilateral contact nonlinearities. Indeed, it is necessary to regularize such nonlinearities using analytic functions. Many analytic functions can be recast into quadratic forms, or defined as solutions of quadratic equations, by introducing (several) auxiliary variables. Hence, the method is not limited to rational polynomials, but also trigonometric, exponential, non-integer power, and logarithm functions can be taken into account, see [START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of nonpolynomial nonlinearities[END_REF] for examples of such recasts. In [109], the application to a piecewise-linear vibro-impact system is demonstrated. The introduction of auxiliary variables results in a (possibly much) larger number of unknowns, which is one of the drawbacks of the method. As discussed in Subsubsection 3.3.5, a crucial aspect of the efficiency of the harmonic balance method in the case of contact problems is the implicit treatment of the strongly nonlinear, and possibly non-smooth contact forces. The quadratic recast requires the explicit treatment of these forces and associated auxiliary variables. To obtain reasonably accurate results, a (much) larger number of harmonics is required to approximate these quantities. This is a possible further downside of the ANM. Lastly, it should be remarked that, to the authors' knowledge, that there is no publication on the application of the ANM to frictional contact problems. Its usefulness in this context thus remains to be seen.

Computation of the branching behavior

At the branching point, one or more new branches emanate from the underlying solution branch, as illustrated in Fig. 13a. Branching points represent one of the main problems encountered during the continuation of the solution in a given parameter interval. The handling of branching points can be divided into detecting that a branching point was passed, localization of the branching point, and switching to the emanating branch(es). Appropriate methods for these tasks are discussed in the following paragraphs. Once a first solution point on an emanating branch is found, numerical continuation can be applied to trace the new branch.

Detection of a branching point A branching point is a bifurcation point, and, hence, the stability of the computed periodic motion changes at this point. A bifurcation detection function is commonly defined, which has a root and changes its sign at the bifurcation point. Two different approaches are often used, one of them is based on the Floquet multipliers, and the other on the Jacobian matrix. One detection approach is to monitor the Floquet multipliers. In the following, we focus on the case of a simple bifurcation point, where a single Floquet multiplier leaves/enters the unit disk in the complex plane. Depending on the way, the Floquet multipliers crosses the unit circle, the type of bifurcation can be inferred [160].

For instance, when a single Floquet multiplier crosses the unit circle at +1, this indicates a saddle node bifurcation, which does not require special treatment. Another detection approach is based on the Jacobian matrix. At the bifurcation point, the Jacobian matrix exhibits a rank deficiency of one (simple bifurcation point). A rank loss of ∂R/∂X indicates a turning point, whereas a rank loss of the extended Jacobian ∂R/∂X T ∂R/∂ρ indicates a branching point.

Localization of a branching point

The solution is typically highly sensitive in the neighborhood of a bifurcation point. Thus, it is often useful to refine the solution branch in the relevant parameter range by continuation with a smaller step size. The location of the branching point can be approximated by interpolation based on the bifurcation test function. More specifically, the point is estimated at which this function actually has a root.

A direct approach for the localization of bifurcation points is to formulate and solve a so-called branching system of equations, the solution of which is the bifurcation point. To this end, Eqs. [START_REF] Ewins | Modal Testing: Theory and Practice[END_REF] are augmented by a constraint equation that enforces that the extended Jacobian has a rank defect of unity. This rank defect is associated with a null space spanned by a single vector, which is considered an unknown of the branching system. The resulting set of equations involves first-order derivatives of the residual. If a gradient-based solver is used, second-order derivatives are therefore required, which are generally more difficult to obtain analytically.

Branch switching Once the location of the branching point is known, the emanating branch(es) can be determined. Two approaches are mentioned here, perturbation and eigenanalysis of the Jacobian. The perturbation approach takes advantage of the fact that bifurcation points are often structurally unstable; i. e. , even a slight perturbation of the problem will unfold the bifurcation and give rise to disconnected branches. For instance, a small random vector can be added to Eqs. [START_REF] Ewins | Modal Testing: Theory and Practice[END_REF] to order to perturb the problem. By using the continuation methods described above, either of the branches will be followed, depending on the perturbation vector.

A more direct approach is based on the eigenanalysis of the Jacobian at the bifurcation point [START_REF] Grolet | Vibration Analysis of a Nonlinear System With Cyclic Symmetry[END_REF]134]. The eigenvector(s) associated with zero eigenvalue(s) represent the tangent(s) of the emanating branch(es). For this approach, it is crucial to ensure a precise localization of the bifurcation point, which can, for instance, be accomplished by the solution of the branching system. With these tangents, the above described predictorcorrector scheme can be used to obtain a first point on the emanating branch(es).

Strategies for finding a first solution point

To start the continuation procedure, a first solution point must be determined. Local methods are commonly be used to compute the first solution point. To start the iteration process of the local method, an initial guess is required. The closer this guess is to the actual solution, the more likely is the successful convergence of the solver and the less iterations are required to achieve convergence.

In some situations, the exact solution is known for a certain value parameter ρ. In a frequency response analysis, for instance, the upper or lower frequency bound might be located in the linear regime, i. e. , where the contact is either sticking or separated; hence, an exact solution can be obtained. In the case of flutter-induced limit cycles, the situation is more complicated, since self-excited limit cycles do not exist in the linear(ized) system. Moreover, the static equilibrium solves the equations of motion in this case, even though this solution is, of course, unstable. Therefore, a particularly good initial guess is required. One approach is a preliminary nonlinear modal analysis, see Subsubsection 3.5.4. This method yields the effective modal damping of the system as a function of the vibration level. By considering the negative modal damping induced by flutter, a reasonable approximation of the flutter-induced limit cycle oscillations can be obtained. More generally, the homotopy method can be used. For this method, a scalar parameter, the so-called homotopy parameter, is introduced. When the homotopy parameter is continuously varied, the original problem is 'transformed' to a problem with known solution, and vice-versa. For example, a scaling parameter for the nonlinear forces can be introduced. If this parameter is set to zero, one has a linear problem with known solution; if it set to one, one has the nonlinear problem of which the solution is actually sought. By varying the parameter iteratively from zero to one, and solving the resulting subproblems, one successively obtains a first solution point of the original problem. Of course, the homotopy method can be equipped with a continuation strategy, where the homotopy parameter is considered as parameter ρ.

In [123], the homotopy method was utilized to determine an initial solution for flutter-induced limit cycles.

To this end, an artificial harmonic external forcing was applied, which is scaled with the homotopy parameter. The homotopy parameter is then successively decreased towards zero, while the forced resonant behavior is computed by means of continuation. At the end of this procedure, the solution of the autonomous system is obtained, namely the flutter-induced limit cycle.

Strategies for finding isolated branches

Besides turning and branching points, isolated branches, as illustrated in Fig. 13a, represent one of the main problems encountered during the computation of the solutions in a given parameter interval. Since the isolated branch is not connected to the main solution branch, it cannot be obtained by continuation and local branch switching from the main branch. It should be noted that once a first solution point on the isolated branch is found, continuation can be used to determine the whole branch. Hence, the main problem is to determine a first point on the isolated branch. Methods that are qualified for this task are referred to as global analysis methods [141].

A simple approach is the multi-start strategy. For this strategy, a number of initial guesses is generated, randomly or in a regular fashion. Starting from every initial guess, or a promising subset of them, a local method is employed to compute the corresponding solutions. The idea behind this approach is that if enough initial guesses are considered, the basins of attraction of each possible solution should contain at least one guess. Consequently, all solutions will be found. If the number of unknowns is small, this is an interesting approach. However, the number of guesses required to cover all basins of attraction grows comparatively fast with the number of unknowns, and can lead to prohibitive computational effort.

For the specific class of polynomial nonlinearities, Eq. ( 26) takes the form of a system of multi-variate polynomial equation. By utilizing the so-called Groebner basis, it is possible to simplify this problem to the problem of sequentially finding all roots of univariate polynomial equations [START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and Groebner bases[END_REF]. For the latter problem, robust methods are available. However, the computation time of this method grows exponentially with the number of unknowns. According to Grolet and Thouverez [START_REF] Grolet | Computing multiple periodic solutions of nonlinear vibration problems using the harmonic balance method and Groebner bases[END_REF], it is quite certain that the method becomes practically infeasible for a number of about 30 generalized coordi-nates. In general, the requirement of polynomial nonlinearities makes the applicability to contact nonlinearities questionable.

In general, the concept of homotopy can also be useful to find isolated branches. Often, the detached branch merges with the main branch if a certain parameter of the problem, such as the excitation level or a contact parameter, is varied. By continuing an appropriate solution point with respect to this parameter, one should eventually land on the detached branch. However, homotopy techniques also suffer from the exponential increase of computational effort with the number of unknowns. An advantage of these methods is that they can be ad hoc parallelized, so that their total run-time can be considerably decreased by implementing the method on parallel computers.

Scaling of unknowns

In general, the different unknowns in Eqs. ( 34) may assume numerical values of different order of magnitude. This is very common if different physical quantities are treated as unknown. For example, this is the case if the frequency is considered as an unknown Ω ∈ X, in addition to the generalized coordinates. This can have a crucial influence on the convergence behavior of the numerical solution procedure. A linear scaling can be applied to the unknowns, so that they have approximately matching orders of magnitude. This can greatly improve the convergence rate and overcome divergence problems.

Advanced simulation techniques

The methods discussed in the previous subsections represent the general framework for the vibration prediction of bladed disks coupled by mechanical joints. Several interesting variants of these general methods have been developed for specific purposes, some of which are addressed in this subsection. In Subsubsection 3.5.1, the direct computation of forced resonances is discussed. In Subsubsection 3.5.2 and Subsubsection 3.5.3, methods are presented that permit the consideration of parameter uncertainty and time-dependence of parameters, respectively. Finally, the concept of nonlinear modes is discussed in Subsubsection 3.5.4, which can be useful for the characterization of the nonlinear vibration behavior and provides a means of model order reduction.

Direct computation of resonances

In the case of forced response analyses, the maximum vibration response within a certain excitation frequency range (resonance) is typically of foremost importance in the design process. Specifically, it is often relevant to determine how the resonant response varies with certain parameters. A straight-forward approach for the parametric analysis of resonances is to compute the frequency response for a discrete set of parameter points, and then extract the resonances from the results. However, this strategy leads to a relatively large computational overhead, since most of the determined frequency response points are not further relevant. Furthermore, the resonant behavior may exhibit a strongly nonlinear dependence on the parameters, so that regions of comparatively low and high sensitivity may exist. As these parameter regions are not a priori known, the appropriate sampling of the parameter points is not trivial and a manual, iterative refinement procedure might be necessary to obtain satisfying results.

A more sophisticated approach is the direct parametric analysis of resonances. To this end, the set of equations governing the forced vibrations (Eqs. ( 27)-( 28) is augmented by a constraint that ensures resonant behavior, the so-called resonance criterion. Moreover, the excitation frequency is no longer treated as a parameter, but it becomes an unknown (resonance frequency). As resonance criterion, the horizontal tangent condition [120] can be used, that enforces a vanishing derivative of a properly defined amplitude with respect to the excitation frequency, which is the necessary condition for a maximum response. The resulting set of equations can be solved with the methods described in the previous subsection. Compared to the crude computation of frequency responses, the direct method only computes the relevant resonance points, which can lead to considerable computational savings. Furthermore, by means of continuation, the discretization with regard to the parameters can be automatically adjusted, such that regions of low and high sensitivity are not a notable problem for this method. Due to its local character, however, the method is not capable of detecting additional maxima in the considered frequency range. This is an important limitation, since the emergence of secondary maxima is a typical phenomenon in the presence of nonlinear modal interactions. Consequently, it cannot be ensured that the traced point is actually the global maximum in the considered frequency range 22 . Therefore, the method is limited to the case of isolated resonances [START_REF] Krack | Multiharmonic Analysis and Design of Shroud Friction Joints of Bladed Disks Subject to Microslip[END_REF].

Consideration of parameter uncertainty

Parameter uncertainty refers to the limited knowledge of the exact value of the parameters used in the physical model. Uncertain parameters occur in every domain of the model. The properties of the underlying structure can be uncertain, such as its local stiffness, inertia and damping properties. The aeroelastic properties can be uncertain, such as the external force level and distribution, and the aero-mechanical coupling. Finally, the contact parameters can be uncertain, such as the friction coefficient, and the initial normal pressure and clearance distribution.

To account for parameter uncertainty, a first step is usually a sensitivity analysis. This step is useful to estimate the overall influence of the uncertain parameters on the relevant output quantities (e. g. resonant response levels, flutter stability boundaries, etc. ). Moreover, this is helpful to determine which of the parameters are most relevant, and which can be neglected in possible proceeding investigations. To determine the sensitivities of output quantities in the neighborhood of a nominal set of parameters, their derivative with respect to the parameter can be analyzed. The sensitivity of the solution can be obtained analytically, if one has nothing better to do, by means of implicit derivation of Eqs. ( 27)-( 28) [121,[START_REF] Krack | A High-Order Harmonic Balance Method for Systems With Distinct States[END_REF]. The chain rule has to be applied accordingly, as the output quantities are a function of the solution. By using these sensitivities, the output quantities can also expanded in a Taylor series, typically of first or second order, depending on the order of available sensitivities. In regimes of strongly nonlinear parameter dependence, these low-order expansions are only useful in a very small range. To determine the parameter dependences in larger ranges, re-computations for a discrete set of parameter points are generally required to obtain reliable results. It is common to describe uncertain character of parameter values in terms of probability density functions. In this case, probabilistic methods can be used to determine the stochastic characteristics of the output quantities (expectation values, standard deviations and higher-order moments). The Polynomial Chaos Expansion represents a suitable framework for the uncertainty quantification. The computation of the expansion coefficients involves the evaluation of multi-dimensional integrals. These integrals can be approximated by means of sampling (or collocation) techniques such as Monte Carlo or sparse grid methods [128]. The advantage of sampling techniques is that they are non-intrusive, in the sense that they rely only on deterministic evaluations such that the vibration prediction method does not have to be modified [119,117]. Intrusive approaches have also been developed, which incorporate the probabilistic integration into the simulation method, see [START_REF] Didier | Nonlinear vibrations of a mechanical system with non-regular nonlinearities and uncertainties[END_REF] for an application to the harmonic balance method.

Consideration of the time-dependence of parameters

So far, we assumed the time invariance of the equations of motion. This is a limitation, since this does not capture effects such as the variability of the operating conditions, or the variability of the interior and exterior mechanical properties of the structure due to wear and tear. Fortunately, all these effects take place on time scales that are typically much larger than the time scale associated with the structural vibrations. This can be exploited by separating the time scales using an averaging formalism [111]. The evolution of the slow dynamics (occurring on the longer time scale)can be described by differential equations in the averaged vibrational quantities. The time evolution on the longer time scale can then be analyzed using numerical integration. The averaged vibrational quantities are determined using the vibration prediction methods described above, considering the slowly varying quantities as constant in each time step. For example, Salles et al. [137,138] investigated the effect of fretting wear in the root joints of bladed disks on their vibration behavior. To this end, the wear process was described by Archard's differential law, and computed by means of numerical integration. It was assumed that only the wear depth changes, which is reflected in the local clearances in the contact interface. The vibration behavior was computed using the harmonic balance method.

Nonlinear modes

Similar to its linear counterpart, the concept of nonlinear modes 23 can be utilized to extract the essential vibration characteristics of nonlinear systems, and to develop means of model order reduction. For recent overviews on the topic in general, we refer to [START_REF] Kerschen | Nonlinear Normal Modes, Part I: A Useful Framework for the Structural Dynamicist: Special Issue: Non-linear Structural Dynamics[END_REF]166]. 23 In literature, the term Nonlinear Normal Mode (NNM) is quite common. However, the term 'normal' may lead to the wrong conclusion that nonlinear modes are orthogonal to each other. Apparently this term goes back to Rosenberg [136], who defined nonlinear modes as vibrations in unison, i. e. where all material points cross their equilibrium points and their extremum points simultaneously. For this type of vibration, the motions take place on so-called modal lines in the generalized displacement space which are normal to the surface of maximum potential energy [166]. However, this property is only valid for symmetric conservative systems, whereas non-trivial phase lags among the oscillations of the coordinates may exist in general. Hence the term 'normal' in this context is avoided in this article.

In the field of bladed disks coupled by mechanical joints, the concept is rather novel and far less established than conventional forced response and flutter analyses. In this subsubsection, we briefly introduce the fundamental notions of nonlinear modes and discuss how this concept can be applied for the benefit of analysis and design of friction-damped blade disks. First of all, what are nonlinear modes? Like in the linear case, they describe the autonomous behavior of the system, as opposed to the behavior under operating conditions described by Eq. ( 2). Hence, the excitation terms are not considered in the nonlinear modal analysis, by excluding the external forces and possibly the self-excitation terms describing unstable aeroelastic effects. Nonlinear modes can then be defined as periodic motions of this autonomous system. However, the autonomous vibrations of nonconservative systems, such as friction-damped bladed disks, are generally not periodic, but they decay with time. To account for their time decay, Laxalde and Thouverez [START_REF] Laxalde | Complex Non-Linear Modal Analysis for Mechanical Systems Application to Turbomachinery Bladings With Friction Interfaces[END_REF] proposed to define nonlinear modes as pseudo-periodic motions, socalled complex nonlinear modes. To compute them, the ansatz of the harmonic balance method in Eq. ( 25) is modified by introducing a real part in the exponential function in addition to the imaginary part related to the oscillation. Here, it is assumed that the damping terms are frequency independent, rendering the approach inexact. Also, inaccuracy is introduced by the approximate evaluation of the nonlinear forces based on the underlying periodic forms. The concept of complex nonlinear modes appears reasonable if the strictly autonomous, decaying behavior is of interest. As discussed in Subsection 3.1, however, the dynamic regime of foremost importance is still that of periodic vibrations. Therefore, it was suggested by Krack [START_REF] Krack | Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept[END_REF] to define nonlinear modes of nonconservative systems as periodic motions, induced by a (usually) negative damping term that is just large enough to compensate the inherent dissipation of the system (periodic motion concept). The results obtained for both methods are similar if damping is comparatively weak [START_REF] Krack | Nonlinear modal analysis of nonconservative systems: Extension of the periodic motion concept[END_REF]. In general, the definition as periodic, undamped motions is better suited to reflect the vibration behavior in the presence of a sustained excitation source (external forcing or negative damping). Moreover, this concept permits the use of standard methods for the computation and stability analysis of periodic motions, including the conventional harmonic balance method and the shooting method. This is an important benefit compared to complex nonlinear modes, which can only be computed by means of the generalized Fourier method. Using either of the above mentioned definitions, the nonlinear modal analysis is carried out. What are the results and how can they be utilized?

(a) Modal properties as a function of the vibration level : For the considered mode, the natural frequency, the modal damping ratio, and the modal deflection shape are obtained as a function of the vibration level. In the case of flutter, possible limit cycles can be determined as the points where the total damping equals zero. In the case of harmonic forcing, the modal properties characterize the resonant behavior of the system according to the deformationat-resonance hypothesis. In particular, the modal properties determine how the resonance frequency shifts, how effectively the system is damped, and if the vibration form changes or even localizes due to nonlinear effects. It should be emphasized that these results are central to the design of bladed disks. In fact, the primary purpose of conventional forced response analyses is often the resonance frequency shift and to determine an effective damping measure. (b) Identification of nonlinear modal interactions: Nonlinear modal interactions refer to the energy exchange between two or more modes of vibration, which is caused by nonlinear effects. Thus, a mode that is not directly forced by an external source, may be excited internally by the action of nonlinear forces. In the presence of modal interactions, the dynamic behavior is usually highly sensitive and can considerably deviate from the behavior without modal interactions. In the near-resonant steadystate frequency response, modal interactions often lead to the emergence of secondary maxima, and appear to be associated with the occurrence of isolated branches [START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF]. In the context of friction damping, numerical and experimental evidence of nonlinear modal interactions was reported e. g. in [170, [START_REF] Berthillier | Blades Forced Response Analysis with Friction Dampers[END_REF][START_REF] Chen | Periodic Forced Response of Structures Having Three-Dimensional Frictional Constraints[END_REF][START_REF] Chen | Periodic Response of Blades Having Three-Dimensional Nonlinear Shroud Constraints[END_REF][START_REF] Krack | A Method for Nonlinear Modal Analysis and Synthesis: Application to Harmonically Forced and Self-Excited Mechanical Systems[END_REF]. Modal interactions can introduce a considerable higher harmonic vibration content, which turns out to have a largely detrimental effect on the effective damping [START_REF] Krack | On the Efficacy of Friction Damping in the Presence of Nonlinear Modal Interactions[END_REF]. Therefore, nonlinear modal interactions are relevant to the design of bladed disks coupled by mechanical joints. The nonlinear modal analysis is well-suited to detect and characterize such interactions with regard to the participating modes and frequencies. Moreover, it can be inferred from the results in which regimes the vibration behavior can be represented by an isolated nonlinear mode. (c) Modal reduction: As in the linear case, the nonlinear modes can be utilized as a basis for the approximation of the (non-autonomous) vibration behavior under operating conditions. Since the superposition principle is not valid in nonlinear systems, the approaches are largely limited to single-modal approximations. The range of validity of such approximations is the regime in which the vibration behavior can be accurately represented by an isolated nonlinear mode (absence of nonlinear modal interactions). In this case, the problem can be reduced to a single-degree-of-freedom oscillator, whose properties are those of the considered nonlinear mode. Note that since the vibrational deflection shape depends on the vibration level, the contributions of all considered linear modes are accounted for in a single nonlinear mode. The computational effort required for the evaluation of this reduced problem is almost negligible compared to the initial problem. This makes the approach interesting for comprehensive parameter studies, design optimizations and uncertainty analysis. For applications to forced response and flutter-induced limit cycle analyses of bladed disks coupled by shroud joints and underplatform dampers, see [START_REF] Krack | A Method for Nonlinear Modal Analysis and Synthesis: Application to Harmonically Forced and Self-Excited Mechanical Systems[END_REF][START_REF] Krack | Reliability Optimization of Friction-Damped Systems Using Nonlinear Modes[END_REF]. In the case of nonlinear modal interactions, it is generally possible to consider multiple nonlinear modes and to utilize a nonlinear projection method. To the authors' knowledge, however, such an approach has not been applied to jointed structures.

Summary and future work

In Fig. 16, an overview of the overall vibration prediction procedure is illustrated. The different steps are connected via signal flows indicating the exchanged physical and computational quantities.

The vibration prediction methods presented in this work represent only the current state of the art of in the field of bladed disks coupled by friction joints. This field is and will most certainly still be a lively field of research in the future, since the success and structural mechanical integrity of next generation's turbomachines depends crucially on the quality of the vibration prediction. Some of the possible topics of future research are listed below.

-Experimental validations, in particular with regard to the contact modeling approaches. -Fluid-structure interactions are currently described in a comparatively simplified way. It remains to be seen, if and how nonlinear aeroelastic effects can be accounted for in the vibration prediction of nonlinearly coupled bladed disks. The development of fully-coupled approaches that take into account both mechanical and aerodynamical nonlinearities, whether -Interaction between vibrations and damaging mechanisms are only rarely addressed. An important example is the interplay between the nonlinear vibrations of the bladed disk and fretting in the joints. An interesting topic could be the initiation and propagation of cracks of vibrating bladed disks with friction joints. -Uncertainty of various system parameters, including those associated with the contact formulation. -Variability of operating conditions is an often neglected aspect that could prove to be relevant in the design of bladed disks. -Multi-stage effects are known for their potentially considerable effect on the vibration behavior; however, they are so far disregarded in nonlinear vibration analyses. -More robust strategies for the continuation, stability and bifurcation analysis are certainly desirable, since the possible branching behavior of the solution is commonly ignored and the stability is only rarely analyzed.

Compliance with Ethical Standards

Conflict of Interest: The authors declare that they have no conflict of interest.

A The traveling wave coordinate system

In this appendix, we define the traveling wave coordinate system and illustrate the notion of traveling waves. The traveling wave coordinates are related to the physical coordinates by the (inverse) discrete Fourier transform. The transformation can be applied to any physical quantity (displacement, force, etc. ). For the displacement, this transform reads

tw u j (t) = 1 √ n s n s -1 n=0 e -i 2πjn ns (n) fe u(t) , j ∈ [0, n s -1] , ( 39
) (n) fe u(t) = 1 √ n s n s -1 j=0 e i 2πjn ns tw u j (t) , n ∈ [0, n s -1] , (40)
tw u(t) = W H n s ⊗ I n fe,s fe u(t) ⇔ fe u(t) = W n s ⊗ I n fe,s tw u(t) . (41)
Herein, tw u j denotes the displacement vector associated with the (spatial) wave number24 j, and H denotes the Hermitian transpose. A congruent discretization and ordering is here assumed for each sector, such that the physical displacement vector

(n) fe u comprises the same number of degrees of freedom, n fe,n = n fe,s for each sector n ∈ [0, n s -1]. The vector tw u j has the same number of coordinates as the physical displacement vector (n) fe u. The more compact notation in Eq. ( 41) involves the Fourier matrix W n s , defined in Eq. [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF].

Remark At this point, it is important to note that both fe u(t) and tw u(t) are complex-valued quantities in general. The complex arithmetic is very convenient for the mathematical derivations. Of course, eventually we are interested only in the physical part, that is, the real component.

To illustrate the traveling wave character of the coordinate system, regard the k-th wave component, tw q k , of a physical quantity q (displacement, force, etc. ), and consider the case of an oscillation with tw q k (t) = √ n s Qe iφ(t) .

Herein, Q ∈ C n fe,s is a complex-valued amplitude vector and φ(t) is the phase, which is assumed to be strictly monotonous in time with φ(t) > 0 in the considered time span, see Fig. 17b.

. . . As a consequence, the relation between phase and time is bijective. All other wave components are assumed to be zero, tw q j = 0 ∀j ∈ [0, n s -1] \k. Taking into account the transform defined in Eqs. ( 39)-( 41), the response of sector n reads

(n) q (φ (t)) = 1 √ n s e i 2πkn
ns tw q k (φ (t))

= Qe i φ(t)+ 2πkn ns = (0) q (φ(t) + nθ k ) .

Herein, the abbreviation θ k is used, with

θ k = 2π n s k . ( 44
)
θ k is the phase lag between neighboring sectors and is hence often referred to as inter-sector/-blade phase angle (IBPA). Since the relation between phase and time is considered bijective, Eq. ( 43) defines a unique time lag ∆t with (n) q(t) = (0) q(t + ∆t). Moreover, q(t) is spatially periodic; i. e. , (n+n s ) q(t) = (n) q(t), which can be easily verified from Eq. ( 43). Hence, the spatiotemporal form of the quantity p(t) can be identified as a traveling wave, discrete in space and a continuous in time, as illustrated in Fig. 17a. Since φ is allowed to be time-dependent, the time lag ∆t generally varies with time t and sector number n. This means that the wave does not have to propagate with constant speed. However, the case of constant wave speed is of particular interest, and it is coincides with a harmonic oscillation of constant angular frequency Ω, φ(t) = Ω. In this case, the constant wave speed is Ω/k (in rad/s).

It should be noted that Eq. ( 43) defines a strict backward traveling wave. However, the traveling wave nature is seen to alias relative to n s , depending on the wave number k. The ranges of k and θ k that correspond to apparent forward and backward traveling waves (FTW and BTW, respectively), and the special cases of standing waves (SW), are given in Tab. In this sense, the columns of the Fourier matrix W n s in Eq. ( 41) can be interpreted as discrete unit traveling waves, such that tw u k corresponds to a wave form with wave number k and IBPA θ k .

Remark It should be emphasized that wave numbers and nodal diameter numbers are only illustrative expressions for strictly mathematical concepts. The physical number of waves or nodal diameters can deviate from the mathematical one.

Consider the example of a rotationally symmetric disk. We can divide the disk into a finite number of n s sectors. The highest possible wave number is then bounded by n s -1 in accordance with our definition. But, of course, the disk can carry an infinite number of waves. The higher wave forms are generally not lost by the dissection into a finite number of sectors, but represent higher modes of vibration for a specific mathematical wave number.

In the above considerations, the complex-valued amplitude vector Q is assumed to be constant in time. If Q(t) depends on time, the strict relation (n) q(t) = (0) q(t + ∆t) between time and sector number is no longer satisfied. However, it still holds that (n) q(t) = (0) q(t)e inθ k ; that is, there is a constant phase lag between the individual sectors for a given k. The corresponding motion for time-dependent Q(t) can thus be interpreted as a pseudo-traveling wave. This notion can be useful to describe vibration phenomena during run-up or run-down of a rotating machine. In this case, the oscillation frequency and the amplitudes vary with time, but excitation and vibration response might still exhibit a characteristic traveling wave form. For instance, in the case of constant acceleration, the phase φ would be defined as φ(t) = α.

The discrete Fourier matrix W n s for a number of n s sectors (or samples in general) is defined as,

W n s = [W ab ] = w 1 • • • w n s , with W ab = w (a-1)(b-1) n s . ( 47
)
where w n s = e

i 2π
ns is the n s -th root of unity.

B Traveling wave structural matrices

In this appendix, we describe how the structural matrices in traveling wave coordinates tw K k , tw D k and tw M k can be obtained from the structural matrices of a reference sector. Suppose that a finite element model of the reference sector is given. The sector spans an angular region of 2π/n s , and can be divided into an inner volume, and left and right boundary 25 . Typically, the finite element model spans the whole sector, including left and right boundaries, making the description of the reference sector somewhat redundant. The sector's displacement vector in physical coordinates can be permuted and partitioned as u T l u i u T r T , where u l and u r are degrees of freedom associated with left and right boundaries and u i are the inner degrees of freedom. u l and u r shall have the dimension n b and u i shall have the dimension n i . Because of the redundancy, the total number of degrees of freedom of the sector (without constraints on left and right boundaries) is 2n b + n i and exceeds the number n fe,s = n fe /n s of degrees of freedom per sector (of the full model) by the number of degrees of freedom of one boundary, n b . The accordingly ordered structural matrices (in physical coordinates) without any constraints on left and right boundaries are denoted

where A ll , A rr and A ii account for the coupling within each boundary and inner volume, and the matrices A li and A ri account for the coupling between boundaries and inner volume. It is here assumed that left and right boundaries are disjunct, so that no coupling exists between them. The matrices tw A k for each IBPA θ k can be obtained as [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling[END_REF],

tw A k = P H k (0) fe AP k , A ∈ {K, D, M } , (49)
with the matrix P k

P k =   I n b ⊗ e iθ k 0 0 I n i I n b 0   . ( 50
)
The resulting matrices tw A k then have the proper size n fe,s × n fe,s . It is assumed in the above formulations that left and right 25 Here, 'left' and 'right' is meant with respect to the rotation axis (and, consequently, the numbering of the sectors) boundary have matching nodes and the local coordinate systems are accordingly adjusted. If this is not the case, Eq. ( 50) has to be adjusted in a straight-forward manner accounting for the coupling of non-conforming meshes and coordinate transformation. Note that the formal relation tw A k = W H n s ⊗ I n fe,s A W n s ⊗ I n fe,s could also be utilized to obtain these matrices. However, the method described in this appendix is much more efficient, since it involves only a single sector.

C Transformation to relative coordinates at the contact interface

In this appendix, we describe how the transformation to relative coordinates discussed in Subsubsection 2.3.1 is applied. To this end, the global displacement vector u is expressed in terms of relative coordinates c g at the contact interface and remaining coordinates u rem ,

u = B T + N B T c g u rem = Lu (b) . (51)
Herein, + denotes the pseudo-inverse, and N A = Null (A) denotes the nullspace of matrix A26 .

In general, the transformation matrix L could be expensive to compute, since it involves the computation of the pseudoinverse and the null space of the comparatively large matrix B T (dimension n d × 3n c , where n d could be n fe,s or n fe ). The computational cost can be considerably reduced by taking advantage of the local nature of the problem, i. e. , by computing the sub-matrices of L separately for each interface. For convenience, the coordinate vector u is rearranged in such a way that the nodal DOFs associated with a particular interface are grouped together. The matrix B then takes the form B = bdiag{B 1 , . . . , B n if } where B n is the local coupling matrix of interface n and n if is the number of interfaces. The n int interior DOFs are not associated with any of the interfaces and form the rear part of the re-ordered vector u, where typically n int 3n c . The null space associated to the interior DOFs is trivial and does not have to be computed explicitly. The matrix L can then be assembled as

L =       B T 1 + • • • 0 N B T 1
0 • • • B T n if + 0 • • • N B T n if 0 0 • • • 0 0 • • • 0 I n int       . ( 52
)
The coordinate transform is applied by substituting Eq. ( 51) into the equations of motion and left-multiplication by L T . It should be noted that the matrix L defined in Eq. ( 52) has full rank, and, thus, Eq. ( 51) defines an invertible coordinate transform.

Resulting structure of the contact force vector The structure of the contact force vector depends on the variant pursued for the definition of the coupling DOFs. If the conventional variant (a) is used, the contact force vector f c takes the form,

f (a) c [u] = Bλ B T u , (53)
where λ represents the actual contact force law formulated in terms of the contact gaps (and/or velocities). The gaps are determined by means of the transform B T u, every time when the contact force vector λ is evaluated, and a multiplication by B is necessary to determine the force vector f c acting on the global displacement vector u. In the case of variant (b), this transformation is applied, once and for all, during the dynamic substructuring procedure. The contact force vector thus becomes

f (b) c u (b) = L T f (a) c [Lu] = L T Bλ B T Lu (b) = I 3n c 0 λ I 3n c 0 u (b) = λ [ c g] 0 . ( 54
)
Owing to the preliminary coordinate transformation, the global contact vector depends and acts on only the first 3n c components of the coordinate vector. Hence, no transformation is necessary during the nonlinear dynamic analysis.

D Craig-Bampton and MacNeal-Rubin method

In this appendix, explicit expressions are given for the matrix T of component modes for the well-known CB and MR methods, see e. g. [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF]. Consider an initial model with a hermitian, positive-definite stiffness matrix K = K H > 0 and a hermitian, positive-definite mass matrix. The associated vector of coordinates u is of the form u = u T ret u T del T where u ret and u del denote the coordinates to be retained in the reduced model and those that are (deleted and) only approximated in terms of generalized coordinates, respectively. In the case of the CB method, the reduction basis is spanned by constraint modes and a set of fixed interface normal modes,

T cb = I n ret 0 n ret ×n del Ψ Φ fixed . ( 55
)
The first hyper-column represents the constraint modes, which are static deformation shapes for a unit displacement applied to one of the coupling DOFs, while the remaining coupling DOFs are kept fixed,

K del,del Ψ = -K del,ret , (56)
where K del,ret refers to the restriction of K to the rows associated with u del and the columns associated with u ret and so on. The fixed interface normal modes, assembled in the matrix Φ fixed in Eq. ( 55), are obtained from modal analysis of the system with all coupling DOFs fixed,

Φ fixed = φ fixed 1 • • • φ fixed n mod , K del,del -ω fixed
It should be emphasized that the coupling DOFs are either the nodal or the relative DOFs at the interface, as explained in Subsubsection 2.3.1. The MR method is the complement of the CB method with free interface normal modes. In the case of the MR method, the reduction basis is, thus, spanned by the residual attachment modes and a set of free interface normal mode shapes Φ free , T mr = I n ret 0 n ret ×n del ∆ Φ free del -∆Φ free ret .

The residual attachment modes essentially represent the static deformation shapes for a unit force applied to one of the coupling DOFs, while the remaining DOFs are not loaded.

∆ = R del R -1 ret , R ret R del = K -1 ret,ret K -1 del,ret
Herein, R denotes the residual static flexibility associated with loading of the retained coordinates. The presence of rigid body modes requires special attention [START_REF] Craig | Coupling of Substructures for Dynamic Analysis: An Overview[END_REF]; however, this case is not further discussed in this work. The free interface normal modes are defined as

Φ free = φ free 1 • • • φ free n mod , K -ω fixed j 2 M φ free j = 0 , ω free 1 ≤ . . . ≤ ω free n mod . (60)

E Exact condensation procedure

In this appendix, an exact procedure is presented for the condensation of the harmonic balance equations, which takes advantage of the sparsity of the nonlinear terms. To this end, it is convenient to arrange the equations of motion in such a manner that the nonlinear force and generalized coordinates vectors have the form

f c [u] = λ [ c g] 0 , u = c g u rem . (61)
Herein, f c and u have the dimension n d , whereas λ and c g have the dimension 3n c . We refer to c g as nonlinear coordinates, and to u rem as linear coordinates, since for given c g(t) a linear ODE governs u rem (t). The vector of nonlinear forces is considered as sparse, if 3n c n d . This sparsity is inherited by the harmonics Λ and the associated gradients. The extent of this sparsity depends on the choice of the generalized coordinates. If the physical coordinates c g, that describe the (relative) interface motions, are not retained, the sparsity is generally lost. In the simplest case, c g represent the local relative deformation at the contact points. Taking advantage of this sparsity during the numerical solution process is a common procedure in conjunction with harmonic balance, see e. g. [START_REF] Kim | Periodic Response of Multi-Disk Rotors With Bearing Clearances[END_REF][START_REF] Groll | The Harmonic Balance Method With Arc-Length Continuation in Rotor/Stator Contact Problems[END_REF][START_REF] Berthillier | Blades Forced Response Analysis with Friction Dampers[END_REF][START_REF] Chen | Periodic Response of Blades Having Three-Dimensional Nonlinear Shroud Constraints[END_REF]. To this end, one condenses the set of n d nonlinear algebraic equations for each harmonic to a set of 3n c equations. This can significantly reduce the number of explicit unknowns and thus reduce the computational effort required for the iterative solution process. The procedure is exemplified for the balance of generalized displacement given in Eq. ( 28), but a fully analogous procedure is available for the balance of generalized forces given in Eq. [START_REF] Coudeyras | A New Treatment for Predicting the Self-Excited Vibrations of Nonlinear Systems with Frictional Interfaces: The Constrained Harmonic Balance Method, with application to disc brake squeal[END_REF], see e. g. [132]. To this end, Eq. ( 28) is split into the individual harmonics, U n + H n (Ω) F c,n (U 0 , . . . , U H ) = U ae,n ∀n ∈ H .

The matrices H k are also partitioned as in Eq. ( 61),

H k = H nl,nl k H nl,l k H l,nl k H l,l k . ( 63
)
With this, the first hyper-row of Eq. ( 62) reads

G k + H nl,nl k (Ω)Λ k (G 0 , . . . , G n h ) = 0 ∀ n ∈ H , (64)
where H nl,nl k is a portion of the matrix H k . Eq. ( 64) only depends on the harmonic components G 0 , . . . , G n h of the nonlinear coordinates, but not on U rem,k , the harmonic components of the linear coordinates. It is thus sufficient to solve Eq. ( 64), which is of much smaller dimension than Eq. ( 62) if 3n c n d . Upon solution of Eq. ( 64) for G k , the remaining portion of the generalized coordinates can be determined using the explicit formulation U rem,k = -H l,nl k Λ k . It should be noted that this dynamic condensation procedure is mathematically exact, so that it does not suffer from poor accuracy like, e. g. , the static (Guyan) condensation procedure. Note that the dynamic compliance matrix is defined as the inverse of the dynamic stiffness matrix. Computing H n by matrix inversion, however, would be time consuming. This is particularly true since H n depends on Ω and, thus, typically has to be re-computed in every iteration. As long as the dynamic stiffness matrix can be expressed as a polynomial in Ω with constant coefficient matrices, the matrix inversion can be replaced by a small number of matrix multiplications and the trivial inversion of a diagonal matrix, see e. g. [132,[START_REF] Krack | A Method for Nonlinear Modal Analysis and Synthesis: Application to Harmonically Forced and Self-Excited Mechanical Systems[END_REF].

Keywordsf 1 Fig. 1 :

 11 Fig. 1: Examples for High-Cycle-Fatigue failures of bladed disks: (a)-(b) debris of the first stage aero engine compressor rotor and its blades [100], (c) fracture of an aero engine low pressure turbine rotor blade[START_REF] Bhaumik | Failure of a low pressure turbine rotor blade of an aeroengine[END_REF]

Fig. 2 :

 2 Fig. 2: The schematic Campbell diagram illustrates natural frequencies, frequencies of synchronous excitation (integer multiples mΩ rot , depicted as so-called load lines), and representative vibration regimes

Fig. 3 :

 3 Fig. 3: Common types of friction joints [143]: (a) roots joints, (b) tip shrouds, (c) underplatform dampers, (d) damper wires, (e) damper pins

Fig. 7 :

 7 Fig. 7: Common contact laws: (a) unilateral laws for the normal contact, (b) friction laws for the tangential plane illustrated by their steady-state hysteresis for one-dimensional harmonic displacement and constant normal load

  (a) nodal coordinates: The coupling coordinates describe the (absolute) deformation of both surfaces meeting at the contact interface. (b) relative coordinates: The coupling coordinates only describe the relative deformation (contact gap) of the contact interface, usually in the local contact coordinate system.

Fig. 8 :

 8 Fig. 8: Classification of coordinates for structures coupled by nonlinear interfaces: (a) coupling based on nodal coordinates, (b) coupling based on relative coordinates

Fig. 10 :

 10 Fig. 10: Categories of excitation mechanisms and resulting vibration regimes: (a) excitation mechanisms, (b) vibrations in the presence of periodic loading

Fig. 13 :

 13 Fig. 13: Important phenomena of a nonlinear dynamical system: (a) conceptual bifurcation diagram, (b) Poincaré section of possible attractors

(a)

  elastic formulation: Friction is modeled by the elastic Coulomb law, normal contact is modeled by the unilateral-elastic contact law. (b) rigid formulation: Friction is modeled by the setvalued Coulomb law, normal contact is modeled by the set-valued unilateral contact law.

Fig. 16 :

 16 Fig. 16: Schematic illustration of the signal flow of the vibration prediction

Fig. 17 :

 17 Fig. 17: Illustration of the spatiotemporal nature of a traveling wave: (a) q as a function of the sector number and the temporal phase similar to [113], (b) temporal phase as a bijective function of time

j 2 M 1 ≤

 21 del,del φ fixed j . . . ≤ ω fixed n mod .

  1. The contact modeling in the mechanical joints is addressed in Subsection 2.2. The reduced description of the system in terms of component modes follows in Subsection 2.3. Finally, aeroelastic effects are addressed in Subsection 2.4.

	2.1 Description of the rotating structure
	2.1.1 Governing equations in the rotating frame of
	reference

  3. Herein, s - n s and s + n s depend on n s , with

	s -n s =	n s 2 -1 n s even n s -1 2 n s odd	(45)
	s + n s =	n s 2 + 1 n s even n s +1 2 n s odd	(46)

Table 3 :

 3 Correspondence between wave number k, IBPA θ k , nature of the apparent wave form and the short name (FTW: forward traveling wave, BTW: backward traveling wave, SW: standing wave, ND: nodal diameter)

	k	θ k	wave form name
	{0}	{0}	SW	TW 0
	1, s -n s 2 } for even n s { n s s + n s , n s -1	]0, π[ {π} ]π, 2π[	BTW SW FTW	TW -k TW ns 2 TW + n s -k

Note that the actual values of(n) fe f c and (n) fe f a still depend on the independent displacement vectors (n) fe u and are, of course, generally not equal to each other.

[START_REF] Batailly | Study of Component Mode Synthesis Methods in a Rotor-Stator Interaction Case[END_REF] Hence, the part fe K related to elastic forces is sometimes referred to as 'tangent' stiffness matrix

The actual Mortar method does not only involve the discretization by means of contact segments, but is commonly associated with an augmented Lagrangian formulation of the contact laws. Here, 'Mortar-like' only refers to the discretization using contact segments.

u is a vector of generalized displacements of dimension n d × 1. In particular, it can stand for fe u or tw u.

This coordinate transform may also involve reduced descriptions of subdomains in terms of component modes, see Subsection 2.3, in which case u refers to generalized coordinates instead of nodal coordinates. This also includes the special case where one or both of the contacting surfaces belongs to a rigid body, where B takes into account rigid body kinematics.

The actual reduction depends on the discretizaton. In accordance with the notation introduced in the previous subsection, for the variant (a), the number of coupling coordinates equals the number of rows of B containing non-zero elements. For variant (b), the number of coupling coordinates is 3n c , that is, the dimension of the vector of contact gaps c g.

For more information and an illustration of such traveling waves, see A.

Note that one approach to formally derive Eq. (17) is the principle of virtual work.

Since this matrix is formulated in the modal space, it is often referred to as modal aerodynamical influence matrix, and its entries are referred to as modal aerodynamical influence coefficients.

[START_REF] Castanier | Consideration of the Benefits of Intentional Blade Mistuning for the Forced Response of Turbomachinery Rotors[END_REF] Note that for assessing the aeroelastic stability, an alternative to the energy method is to carry out a fluid-structure simulation of the whole annulus, starting from an initial perturbation, and analyzing whether the vibrations grow or decay[165].

As discussed in Subsubsection 2.1.2, this is the case if the sector-to-sector deviations of (geometrical, material and contact) properties are sufficiently small or the inter-sector coupling is sufficiently strong.

In the literature, other widely used names for the method described here are the 'Describing Function method' and the 'Krylov-Bogoliubov-Mitropolsky method'. Moreover, the prefixes 'multi' or 'high-order' are often used for the harmonic balance method in order to clarify the difference to the singleterm variant which only considers the fundamental harmonic.

Hence, this procedure is considered a pure frequencydomain method, since there is no need to switch to the time domain, in contrast to the alternating frequency-time scheme presented in the following paragraph.

[START_REF] Charleux | Threedimensional Multiharmonic Analysis of Contact and Friction in Dovetail Joints[END_REF] cf. discussion in paragraph 'Tangential contact' in Subsubsection 2.2.3

Note that the normal preload is sometimes also referred to as initial normal load, which emphasizes that the actual normal load may change due to vibrations. In fact, even the static component (or average value) of the normal load is influenced by the dynamic contact interactions.

The task of finding such an appropriate guess is addressed in Subsubsection 3.4.4.

In the literature, the name 'Newton-Raphson' method is also commonly used for the method described here.

In fact, in the case of the horizontal tangent condition, it cannot even be ensured that the traced point is a maximum, but a local minimum or a saddle node might be traced 'by mistake'.

In the literature, this number is also referred to as 'harmonic index'. To avoid confusion with (temporal) harmonics in the context of frequency domain methods, this terms is avoided in this work. Moreover, the term 'nodal diameter number' is also common for this number.

Note that B T u = B T B T + B T N B T c g u rem = c g, infull accordance with Eq. (6).