Particle swarm optimization for support vector clustering Separating hyper-plane of unlabeled data

Keywords: Support Vector Machines, clustering, Particle Swarm Optimization

The objective of this work is to design a new method to solve the problem of integrating the Vapnik theory, as regards support vector machines, in the field of clustering data. For this we turned to bio-inspired meta-heuristics. Bio-inspired approaches aim to develop models resolving a class of problems by drawing on patterns of behavior developed in ethology. For instance, the Particle Swarm Optimization (PSO) is one of the latest and widely used methods in this regard. Inspired by this paradigm we propose a new method for clustering. The proposed method PSvmC ensures the best separation of the unlabeled data sets into two groups. It aims specifically to explore the basic principles of SVM and to combine it with the meta-heuristic of particle swarm optimization to resolve the clustering problem. Indeed, it makes a contribution in the field of analysis of multivariate data. Obtained results present groups as homogeneous as possible. Indeed, the intra-class value is more efficient when comparing it to those obtained by Hierarchical clustering, Simple K-means and EM algorithms for different database of benchmark.

I. INTRODUCTION

The Support Vector Machines (SVM) appeared for over fifty years and was proposed by Vapnik [START_REF] Vapnik | The nature of statistical learning theory[END_REF]. They constitute one of the few learning methods that are completely derived from theoretical considerations. Therefore it is mathematically well-founded. Indeed, these large margins separators are based on two key ideas, allowing the processing of nonlinear discrimination problems, and the reformulation of the classification problem as a quadratic optimization problem. Moreover, these methods allow the pooling of the maximum margin concept and the notion of kernel function. SVM are supervised classification techniques based on the principle of structural risk minimization rather than empirical risk.

Thanks to its theoretical basis and to benefit from its used cores, we are very interested to exploit the advantages of SVM for the problem of data clustering. Indeed, the enormous quantities of information gathered from various sophisticated technologies such as computers and satellites, have been born an urgent and insistent need to intelligent tools for the analysis, management and organization of this bulky information.

In another words, the clustering term is associated with the notion of abstraction. It is used to summarize information in groups most generals. Indeed, a clustering aims to bring together a set of observations in homogeneous classes by promoting heterogeneity between these different classes. In this context, the general principle of clustering consists of minimizing the distance between two individuals of the same class and maximizing the distance between two individuals of distinct classes.

Besides in recent years, the living world can offer methods to problem-solving that are inspired by its mechanisms. Indeed, there are varieties of research in the literature showing the ability of bio inspired algorithms to produce good results for complex and large problems, such as clustering. The particle swarm optimization is a new bio-inspired optimization method that presents simplicity of algorithms, a convergence speed and excellent results.

Then originally, the method PSvmC (Particle of SVMs for Clustering) was created specifically to address the problem of using SVM in a different field of classification. This problem, which is the main application of this paper, is a clustering problem.

In this paper, first of all we introduce the various work using SVMs in the clustering field. Then, we present a detailed description of our algorithm, its formal presentation and the various experiments performed.

II. RELATED WORK

The main idea of SVMs is that these methods were developed to solve a classification problem (in two classes) by computing an optimal class separator. The separating hyperplane is calculated by projecting the training set defining the two classes in a higher dimensional space. The classification of each data point will be done according to its position to the separator.

In other terms, the optimal hyper-plane is the hyper-plane whose minimum distance to the training examples is maximal: Maximize margin subject to classifying all points correctly.

SFAX University

Sfax, Tunisia yassine.benayed@gmail.com The Fig. 1 shows the different linear hyper-plane allowing the separation of data and the optimal hyper-plane ensuring the maximization of the margin. Solving the problem of calculating the optimal hyper-plane is done in a completely mathematical and theoretical manner (more details in [START_REF] Chaabouni | Particle swarm optimization for support vector clustering[END_REF]). Thanks to the advantages of SVM methods and their solid theoretical foundation, several researchers have attempted to use it in the clustering field but in a semi-supervised way. Indeed, there are few adaptations for pure unsupervised case.

The most famous clustering algorithm based on SVMs is the one invented by Vapnik and called SVC (support vector clustering). We will detail the SVC principles in the following of this section. Also we will make an overview of the hybridization attempts of SVMs with other clustering methods (hierarchical method or k-means method).

 SVC basing on the approach of SVMs, [START_REF] Ben-Hur | Support Vector Clustering[END_REF] has created a novel clustering method which identifies clusters by the separation of each one with a set of contours. This method proceeds by mapping the data points in a high dimensional feature space with the support of a Gaussian kernel. Afterwards, it searches the smallest sphere that encompasses the data image being in the feature space. Finally, it ensures the mapping back of the found sphere to the data space. That creates a set of contours which enclose the data points. Cluster designs points enclosed in the same contour.

 Support Cluster Machines: This method is proposed by [START_REF] Li | Support cluster machine[END_REF] to effectively deal with large-scale classification problems. The main idea of the SCM is to cluster the training set and then apply the SVM classification method on clusters previously obtained.

 CB-SVM: The classifying large data sets using SVMs with hierarchical clusters (using the algorithm BIRCH) presents a semi supervised classification method proposed by [START_REF] Yu | Classifying Large Data Sets Using SVMs with Hierarchical Clusters[END_REF]. It allows the recursive selection of the clusters centroids as the representatives along the hierarchical clustering tree.

 The cluster-based data selection method [START_REF] Sun | Cluster-based Support Vector Machines in Text-Independent Speaker Identification[END_REF], presents a cluster-based learning methodology to reduce training time and the memory size for SVM. Using the Kmeans clustering technique, a set of cluster is created and then filtered. A set of critical data is so prepared to be classified with the SVM method.

III. DESCRIPTION OF THE PSVMC METHOD PSvmC 'Particle of SVMs for Clustering' is our proposed algorithm. It allows the combination of the PSO meta-heuristic [START_REF] James | Particle swarm optimization[END_REF], with the benefits of SVM: the supervised classification method for clustering purposes.

PSvmC provides grouping of data in two different classes. Indeed PSvmC is an algorithm able to create a separator of data without labellum. To do this, the swarm, turned by the PSvmC algorithm, presents a collection of particles. Their number is well determined. Each of these particles consists of defining an optimal hyper-plane separating the data into two classes. These particles move in the data space, keeping in mind the best separator already found, with a velocity properly adjusted. The moving of all particles of swarm led to discovery a global solution of the problem of finding the best separating of data. Hereafter we will detail the main elements of the PSvmC method.

A. Score function

In our work, the measurement of the performance quality of the solution is made through a fitness function to maximize. We seek to choose an objective function which estimates the group quality that we want to refine and improve. In the literature, the internal validity indices such as the inertial indices are the most objective functions used by the famous clustering algorithms, this is explained by the computing quality and performances offered by these indices.

Therefore, after a study of different criteria that ensure a good clustering [START_REF] Chaabouni | Particle swarm optimization for support vector clustering[END_REF], we choose to use inertial index: Intra-class inertia, because it is the simplest and most effective measure used for calculating the raw similarity between objects of same class. Furthermore, the computation time of the intraclass inertia is reasonable and the noisy data does not produce a calculation problem. This index is in a digital format that we want to minimize:

  intra = 1/n × ΣΣ d 2 (x,g i )



Where n presents the number of elements and g i is the gravity center of the class i.

To receive the benefits of this measure we should made a clear expression of distance between objects. We choose the Manhattan distance because firstly we have sets of multivariate data and secondly this measure is not sensitive to noisy data. Indeed, the Manhattan distance calculates the sum of absolute values of differences between the coordinates of two points:

 d (X, Y)= Σ | x i -y i |  B. Particles
The particle notion originating of nature presents a simulation of bees or birds that use suitably chosen paths in their problem of foraging. A particle is then moving, and it has a velocity. Also, each particle has a small memory, which allows remembering its best performance and position. Also, it W Margin SV SV Ho H has a group of informants, generally called its neighborhood. In our method, we consider that the particles are characterized by four features:  A current separating hyper-plane presents the current position of the particle.

 Current intra-class inertia of data which describes the performance value, also called the objective function at the current position.

 The best hyper-plane and its best intra-class inertia, previously found by the particle itself. This is the simulation of the memory characterizing a particle.

 The best hyper-plane and its best intra-class inertia, discovered by the entire swarm. This is to emphasize the neighborhood topology chosen.

To better explain the concept of particle, we will detail in the following: the position, the objective function, the particle velocity and the rules for updating velocity and mutations submitted to renovate the particle swarm.

C. Particle position

In this work, the particle position presents a hyper-plane formed by an array of support vectors which are generated by the SVM method. The support vectors "SV" are the data samples that maximize the distance or the margin separating the positive samples from the negative ones, or simply separate the two classes of the database. The support vector synthesizes someway the important aspects of the sample set. In fact, only the examples corresponding to support vectors are actually useful in learning. If we knew them a priori, we could make learning without considering of other examples. Briefly, in our program we used model files which characterize hyper-planes, theoretically written as follows,

 h(x)=(x)+b, 

Where  is the weight vector  = ( 1 … n ) and b is a constant. A model file contains a set of support vectors and some parameters such as the constant that specifies the maximum-margin between the data-points of the two classes noted b or rho.

A particle at a certain position is a set of SV at a given iteration. This means that a movement of a particle causes a change of the SV set of the hyper-plane.

D. Particle velocity and its updating rules

Let V t i , be the velocity of the particle i at iteration t. The new velocity of i is calculated as a linear combination of three elements by using confidence coefficients:  The adventurous tendency, i.e. continue under the current speed of the current position "current position ".

 The conservative tendency bringing more or less towards the best position already found "best particle ".

 The tendency of overreaction, pointing towards the best informant position of the swarm "best neighbor ".

The application of the combination of these three elements to the current position gives the new position of the particle. Formally the equation of velocity is written as follows with c 1 , c 2 and c 3 are the confidence coefficients where their sum equals to 1.

V i t+1 = c 1 × V i t + c 2 × (best particle -current position ) + c 3 × (best neighbor -current position ) 

As we have defined the position of particle by the set of support vectors, the velocity at a given iteration will be defined as the combination of a number of support vectors taken from the current position "currentSV", a number of support vectors taken from the best position already found "pBestSV" and a number of support vectors taken from the best position achieved by the swarm "gBestSV".

The simulation of the subtraction (best particle -current position ) (respectively best neighborcurrent position ) returns the support vectors belonging to the pBestSV (respectively gBestSV) and outside of the currentSV. As we are currently managing sets, the subtraction will be written as the set of pBestSV (respectively gBestSV) private of the currentSV set: pBestSV \ currentS V (respectively gBestSV \ currentSV).

The velocity is then written as follows:

V i t+1 = c 1 × currentSV+ c 2 × (pBestSV \ currentSV) + c 3 × (gBest \currentSV)  
In our approach, the correlation coefficients represent the percentage of retained SVs from the result of the previously subtraction operation. The value of these coefficients is not arbitrarily chosen, it follows a definite rule, which we will detail in the following:

-Rules for updating the velocity:

The random selection of the correlation coefficients needs a new parameter describing the maximum velocity of particle which must not be exceeded because generally, a high velocity causes a poor exploration of the search space. Our idea to simplify and minimize the number of parameters to be determined, leads us to create a new rule for determining the correlation coefficients. The latter are determined automatically depending on the intra-class inertia of each particle position. As the constraint:

c 1 c 2 c 3 = 1  
So c 1 , c 2 , c 3 must belong to] 0, 1[, and knowing that inertia is a positive and strictly greater to 1 so its inverse satisfies the constraint and presents a value between 0 and 1. To avoid being trapped in very small values and to normalize the inverse of the inertia, we divided it by a common divisor which presents the sum of the inverses of all intra-class inertia "Ia": That is to say the sum of the inertia inverse of the current position position, the inertia inverse of the best particle position and the inertia inverse of the best neighbor . Divisor =

 

Correlation coefficients are then written as follows:

c 1 = c 2 = c 3 = 1 -(c 1 + c 2 )
; (for satisfying the constraint ( 6)).

Calculating these coefficients ensures that the coefficient value increases when the particle position has much more interesting performance. So, when a position has minimal intra-class inertia, the coefficient value will be maximized and then the particle can take as much information as possible from this position.

-Rules for updating the position After calculating the new velocity, each particle considers this information to his new movement. Indeed, we multiply each correlation coefficient by the total number n tot of support vector of the new position. So the new position consists of n tot ×c 1 support vectors from the current position, n tot ×c 2 support vectors from the best particle position and n tot ×c 3 support vectors from the best neighborhood position. As a certain position is composed by a set of support vectors presenting the crucial and decisive samples of the two-class data, the choice of a K number of support vectors in a given position arise the following question:  How many SV must be chosen from each class?

To answer this question, we used the same rule for determining the correlation coefficients seen above but by using the intra-class inertia on each class. The number of SV selected from each class i, is:

 i = n tot × c k ×  
With k = 1… 3 and I i presents the intra-class inertia of each class i.

The selection of N i support vectors of class i, is done by a completely random manner.

E. Mutation operator

To prevent a premature convergence of a search algorithm, the evolutionary programming paradigm introduces an original concept called mutation. The introduction of a mutation operator provides our method with a better exploration capacity of the search space. Thus, the algorithm can converge to the global optimum. This operator consists in changing the current position randomly with a probability equal to 2%. The mutation operator of the algorithm PSvmC embodies the type of deletion mutation. Indeed, it removes 10% of support vectors of a given position.

F. The swarm

The collective behavior of particles defined above leads to the evolution of the swarm which tries to reach interesting sites. Indeed, the random initialization of the particles positions with uniform distribution provides an excellent exploration of the search space and then more chance to converge to the global optimum.

To further explain this choice, the use of the random criterion allows obtaining good and bad position of particles which ensures the separation of data. In other terms, we allow a great value of intra-class inertia. This variety of quality can have a significant gain in terms of solutions variability. This promotes a convergence towards the correct solution beyond the local optimality.

G. Number of particles needed in the swarm

The quantity of particles assigned to solve a problem depends on two parameters:  The size of the search space.

 The relationship between the computing capacity of the machine and the maximum search time.

We have proposed to determine the particle number automatically depending on the attribute number of the database. Indeed, if N is the attribute number, the size of the swarm is given by:

Nbr particle = (K × N)  
With K represents the determination coefficient of the particles number.

Theoretically, a study is needed to extricate the value of the coefficient K. For this, we performed experiments on a set of two-class databases (see the next sections). For each value of the coefficient K, the program of PSvmC algorithm was run 10 times by returning the value of intra-class inertia. To determine the value of the coefficient K we analyzed the development of the average of the 10 intra-class inertia relative to different values of K: 1, 2, 3… We did get the results shown in Fig. 2, Fig. 3. (See all experiments in [START_REF] Chaabouni | Particle swarm optimization for support vector clustering[END_REF]). As the analysis of these results presents that the values of the intra-class inertia shows good stability starting from the value K = 3 for the different databases of benchmark treated, the number of particles as well adopted is Nbr particle = (3 × N).

IV. FORMAL PRESENTATION OF THE PSVMC ALGORITHM

The global PSvmC process is based on three steps: initialization, update, and displacement movement. These steps are crucial to the process of the PSvmC algorithm. This process is summarized by the Fig. 4. During the updating step, each particle of the swarm memorizes the current position and its intra-class inertia if this latter is smaller than its better inertia even found (pBestposition). Besides, the global position (Gbest-position) is updated by the current position if the constraint number (2) presented in Fig. 2 is satisfied. The moving and the mutation steps are applied as we have mentioned in the above section.

The different steps of our method are described in the Algorithm 1: V.

VI. EXPERIMENTS AND RESULTS

A. Used data sets

The data sets used for testing our approach are readily made available by the University of California at Irvine [UCI]1 and the SVM library LIBSVM2 that contain the most commonly used databases. In our simulations we used a set of different two-class databases, designed for simple supervised classification.

 Pima Indians Diabetes (1): used since 1990, this database presented in eight attributes, helps to diagnose and to determine whether a patient shows any sign of diabetes (500 objects) or a good health (268 objects).

Fig. 1 .

 1 Fig.1. Separation of two sets of points by linear separators

Fig. 2 . 1 treating

 21 Fig.2. Effect of the particle number on the intra-class inertia made on the Liver Disorders database.For the Liver Disorders database, since the value K = 3, the intra-class inertia values are limited in the interval[4676; 4900]. This presents a good value of intra-class inertia. Also in

Fig. 3 .

 3 Fig.3. Variation of the intra-class inertia of the Heart database according to the particle number.

Fig. 4 .

 4 Fig.4. Global processing of the PSvmC algorithm.



  Liver disorders (2): presents a set of individuals who complains of a liver disorder. This database contains 345 instances, each one is presented by seven attributes and there are 145 positive and 200 negative instances.  Heart (3): it is made available by LIBSVM. It contains 270 objects described by 13 attributes and containing two classes (120 instances presents a heart disease and 150 instances shows a healthy cardiovascular system).  Breast cancer (4): this database which contains 683 instances (444 positives and 239 negatives) and characterized by nine attributes, describes if a woman has breast cancer, benign or malignant.  Australian (5): it concerns the requests of Australian credit card. It presents a good mix of continuous, nominal attributes with small number of nominal values, and with a greater number of values. It is described by 14 attributes (383 negative, 307 positive).

Table 1 :

 1 Various tests of inertia calculation.

Table 2 :

 2 Various tests of inertia calculation.

UCI is available in the following web address: http://archive.ics.uci.edu/ml/

LIBSVM

Initialize its Current position (i) position with labels randomly taken and its personal best position best particle Calculate the intra-class inertia to affect it to: Ia currentSV and Ia pBestSV .

B. Evaluation

In this section, we describe the experiments that were performed to illustrate the performance of the proposed algorithm PSvmC. To evaluate our work, we made experiments on different data sets, by varying their sizes and comparing the proposed approach to three other algorithms based on different clustering methods: The hierarchical clustering method, the simple k means clustering algorithm and the probabilistic method "EM".

These different grouping methods are available at the Environment for Knowledge Analysis: WEKA [START_REF] Hall | The WEKA data mining software: an update[END_REF]. To evaluate the different generated clusters, we adopted the intraclass inertia to measure the similarity between the same class data as evaluation criteria. We launched our proposed algorithm ten times for each experiment and we calculate at each time the inertia intra-class average. Thus we calculate the intra-class inertia of the three others clustering methods by using a good setup. In fact, if we put a wrong parameter of class number, the intra-class value of these methods becomes incomparable with our method.

Database

Clustering methods

Hierarchical clustering

Simple K-means EM PSvmC

Avg Min According to the results of the four algorithms presented in the table above (Table3) and tested with the different databases of benchmark, we can deduce various observations. Indeed, we noted that the clusters found by PSvmC have more interesting results than those found by other algorithms. In effect, for the "Pima Indians Diabetes" database, our approach provides a far smaller inertia (24063.95) than those obtained by the other three algorithms. Also, for the database of "Liver Disorders", the table 3 shows that by using PSvmC algorithm the intra-class inertia (3799.20) is also smaller than other methods. For the "Heart" database, we see that the difference between the average of the intra-class inertia and the minimum of the intra-class inertia is low; this shows the stability of our approach throw different executions. Finally, we note that the intra-class inertia found by our approach is more efficient when comparing it to those obtained by Hierarchical clustering, Simple K-means and EM algorithms for the "Breast cancer" database (176.07) and for the "Australian" database.

To explain these obtained results: since we actually define a good clustering as the clustering having the small value of intra class inertia, we turned our algorithm in order to find the smallest value of intra-class inertia. A review of the results of various experiments ensures that we have reached our initial goal.

To conclude, the PSvmC algorithm provides a grouping ensuring a smallest intra-class inertia value than the other existing algorithms. Or, the aim of clustering data is to find small intra-class inertia within the same group to increase the similarities between its elements; we can say that our proposed algorithm achieves better results than the K-means, Hierarchical and EM algorithms.

VII. CONCLUSION

In this paper, we introduced our clustering algorithm called PSvmC allowing the creation of a data separator based on the SVM principles. In effect, the particle swarms are a pertinent source of inspiration for integrating SVM in the field of clustering. Moreover, we could propose a heuristic based on PSO, which aims to solve the problem of knowledge extraction.

We described the various experiments carried out to release the best of our approach and this by the right choice of the parameter values of our algorithm. In addition, several tests were conducted to evaluate the performance of our method using intra-class inertia. A comparison on PSvmC was conducted face the Hierarchical clustering, the Simlpe Kmeans and the EM algorithms.

Thanks to the obtained encouraging results we have the idea of proposing other extensions of the PSvmC algorithm that are a top-down hierarchical clustering and a multiobjective clustering. Also, other prospect is related to the improvement of the presented approach such as the determination of the different adjustment parameters of the SVM methods. Indeed, the best setting may improves over the clustering results. Another perspective is the parallelization of the proposed approach to take maximum advantage of the PSO concept that ensures the parallelization of the optimization search.


Nicolas Salvat

Alain Batailly

Mathias Legrand

Two-dimensional modeling of unilateral contact-induced shaft precessional motions in bladed-disk/casing systems

Keywords: turbomachinery, unilateral contact and friction, nonlinear dynamics, rotor-stator interactions, Lagrange multipliers, whirl, modal coincidence

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In the aeronautical and space industry, reliability and efficiency stand as major concerns due to both economical and environmental reasons. Among other challenges, to gain fuel consumption while increasing power, lowering noise and gas emissions, increasing the durability of components while decreasing weight, have attracted considerable interest. Thus, various strategies have been adopted in distinct fields, such as the use ultralight-weight materials (e.g. titanium alloys and composite materials) [START_REF] Williams | Progress in structural materials for aerospace systems[END_REF], planning cost-efficient flight trajectories [START_REF] Félix Patrón | New methods of optimization of the flight profiles for performance database-modeled aircraft[END_REF] as well as the development of multidisciplinary optimization algorithms [START_REF] Sobieszczanski-Sobieski | Multidisciplinary aerospace design optimization: survey of recent developments[END_REF], to name a few. In the specific area of commercial aircraft propulsion, where thrust is often provided by a gas turbine as the one illustrated in Fig. 1, efficiency enhancement [START_REF] Mattingly | Aircraft engine design[END_REF] is achieved through flow-path and geometry optimization, or compression ratios and combustion temperatures increase. One of the most sensitive parameters to tune in compressor and fan stages is the blade-tip clearance, defined as the effective distance between the rotating airfoil tips and the surrounding stationary casing at nominal operating conditions. In fact, this distance must be kept to a minimum to ensure high compression ratios and engine performance [START_REF] Howard | Performance impact of tip clearance variation on a transonic, low aspect ratio, axial compressor stage[END_REF]. As detailed in [START_REF] Lattime | High-pressure-turbine clearance control systems: current practices and future directions[END_REF], potential structural malfunctions may arise from, for instance, mass imbal-ance or blade loss, impact with external objects, casing distortion due to high thermal gradients, fluid excitations through bearings or aerodynamic loads, among others. Due to the increasingly tight operating clearances, these incidents generally lead to involved shaft/bladed-disk/casing interactions induced by unilateral and frictional contact occurrences between rotating and stationary components. Such events that may threaten the engine structural integrity are commonly referred to as structural rotor/stator interactions [START_REF] Bartha | Dry friction backward whirl of rotors[END_REF][START_REF] Batailly | Numericalexperimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF].

Z o Y o Y o X o
rotor bearing (a) Rotor/bearing [START_REF] Muszynska | Whirl and whip-Rotor/bearing stability problems[END_REF] bladed-disk casing (b) Bladed-disk/casing [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF] Figure 2. Systems subject to structural rotor/stator interactions Within this industrial framework, two types of systems must be distinguished depending on the contact interface involved in the interaction: Rotor/bearing systems The rotor is modeled by a flexible shaft and rigid disks, and supported by journal-bearings (Fig. 2(a)), e.g. magnetic [START_REF] Jiang | Improvement of rotor performance under rubbing conditions through active auxiliary bearings[END_REF] or hydrodynamic [START_REF] Sweet | Non-linear rotor bearing behavior[END_REF] bearings, liquid or gas annular seals [START_REF] Childs | Turbomachinery rotordynamics: phenomena, modeling, and analysis[END_REF], among others. Contacts shall occur between shaft and the inner part of its supporting bearings due to large lateral vibrations. Bladed-disk/casing systems A rigid shaft and flexible blades are commonly considered. Contacts take place between blade-tips and surrounding casings (Fig. 2(b)). The blades distortion and the casing wear are the main concerns [START_REF] Batailly | Numericalexperimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact[END_REF][START_REF] Salvat | Modeling of abradable coating removal in aircraft engines through delay differential equations[END_REF]. Generally speaking, three distinct phenomena arising from rotor/stator interactions are discussed in the literature: (1) modal coincidence [START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF], characterized by an exchange of energy between modes of the rotor and the stator of matching nodal diameter through intermittent contacts that may lead to catastrophic structural failures [START_REF][END_REF], [START_REF] Félix Patrón | New methods of optimization of the flight profiles for performance database-modeled aircraft[END_REF] rubbing, refers to direct contact either between shaft and journal bearings [START_REF] Bartha | Dry friction backward whirl of rotors[END_REF] or blade-tips and surrounding casing [START_REF] Legrand | Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law[END_REF] depending on the system of interest and (3) whirl [START_REF] Muszynska | Whirl and whip-Rotor/bearing stability problems[END_REF], describes precessional orbits of the shaft leading to unstable vibratory motions and is related to rubbing problems once the whirling amplitude is larger than the initial rotor/bearing clearance [START_REF] Bartha | Dry friction backward whirl of rotors[END_REF]. The rotor may undergo forward or backward whirl depending on the whirling direction with respect to the shaft rotational speed.

A majority of the numerical investigations on whirl involve flexible shafts with non-deformable cross-sections and the literature is rather scattered when considering whirling motions in fully flexible bladed-disks/casing systems. One of the only available studies [START_REF] Padovan | Nonlinear Dynamics of Rotor-Blade-Casing Rub Interactions[END_REF] targets the transient response of a bladed Jeffcott rotor system. Three response configurations are reported: (1) no contact, (2) single blade contact and (3) multiple blades contact. Backward whirling motions are observed, maximizing blade displacements and stresses. A recent extension of this model was developed in [START_REF] Ma | A revised model for rubbing between rotating blade and elastic casing[END_REF], proposing a novel explicit formulation of the normal rubbing forces dependent on the system parameters such as the rotational speed, the blade cross-section and disk diameter, which was validated with experimental data from a test rig.

A fully flexible rotor is also developed in [START_REF] Lesaffre | Contact analysis of a flexible bladed-rotor[END_REF] using an energetic approach in the rotating frame and considering frictionless sliding contacts via a linear penalty method. The casing is modeled as an elastic ring and the stability of the two coupled structures is addressed through eigenvalue computations, showing different divergent instabilities and mode couplings. Parent et al. [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF] built on this work, using a mass imbalance at the LP turbine as a source of excitation and highlighted two different contact scenarios: the regular case, where contact occurs in the same side than the shaft displacement and results in backward whirling, and the inverted case, where contact occurs in the opposite side than the shaft displacements which leads to forward whirling. A novel kinematic methodology for the calculation of the blade-to-casing distances was incorporated in [START_REF] Parent | 3D interaction in bladed rotor-to-stator contact[END_REF].

In the present study, attention is paid to the shaft/bladeddisk/casing structural couplings and the emergence of whirling motions in the fan stage of a two-spool commercial aircraft engine (Fig. 1) induced by direct unilateral blade/casing contact occurrences. Two contact initiation mechanisms are investigated: [START_REF] Williams | Progress in structural materials for aerospace systems[END_REF] an initially prescribed casing distortion and (2) a mass imbalance on the bladed-disk. Based on previous work [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF], the interactions between these two flexible structures are studied by means of a time-marching strategy similar to [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF]. The corresponding modeling assumptions-detailed in the following sections-are summarized as: Bladed-disk and casing: flexible 2D in-plane models composed of beams. Shaft flexibility: two linear springs. Frame and bearings: array of linear springs. Damping: modal damping with single damping coefficient. Neglected structural terms: centrifugal, spin softening and gyroscopic effects. Rotational velocity: constant in counter-clockwise direction. Solution method: explicit time-marching. Unilateral contact enforcement: Lagrange multipliers method. Frictional contact forces: linear Coulomb law in sliding.

First, the structural 2D in-plane models of the engine stage and the coupling methodology are presented. In a second section, some of the modal properties of the bladed-disk/casing system are explored, highlighting the influence of the shaft flexibility and the role of the frame/bearings coupling on the overall dynamics of the system. Lastly, the simulation procedure is introduced along with a brief analysis of convergence issues related to the implemented numerical strategy. These are followed by the results and discussion of emerging trends, mainly focusing on the types of operating regimes detected and the localization of the strain energy of the dominant modes involved in the interactions.

Structural models

Bladed-disk and casing models

The fan stage model in this investigation extends the one introduced in [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF] and incorporates the flexibility of the shaft through two linear springs connected to the center node of the disk driving the so-called suspension modes. This modelling innovation targets potential whirl motions as the blade dynamics and the in-plane motions of the shaft are now coupled. The bladed-disk, illustrated in Fig. 3(a), is composed of an assembly of straight beams. Each element has two nodes with three degrees of freedom (DoF) per node and is a combination of simple Bernoulli beams (flexural displacements) and rod elements (longitudinal displacements), allowing to account for both radial and tangential displacements. As in [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF], the model consists of 22 blades, with radius of R bd D 0:5 m and ten elements per blade, as well as a dedicated set of beams representing the disk and allowing the structural coupling between neighbouring blades. The stiffness of the center node linear springs is arbitrarily chosen to k o D 3:4 10 7 N=m, which corresponds to about 500 times the tangential stiffness of a blade and is representative of an industrial fan.

The casing, also depicted in Fig. 3(a), is built using a set of curved beams (40 elements). These elements have two nodes and four DoF per node, leading to a total of 160 DoF. Its radius is defined by the radius of the bladed-disk plus an initial clearance ı arbitrarily prescribed: R c D R bd C ı. The considered material is steel for both structures, the corre-sponding geometry parameters are listed in Tab. 1 and the equations of motion governing the dynamics of the system are:

Ä M c 0 0 M bd  R x c R x bd à C Ä D c 0 0 D bd  P x c P x bd à C ÄÄ K c 0 0 K bd C K  x c x bd à D  f c ext f bd ext à C Ä B c B bd f cn (1)
that is, in a compact form:

MR x C DP x C Kx D f ext C Bf cn ( 2
)
where x is the n-dimensional displacement vector, M is the mass matrix and D is the damping matrix computed from the modal domain 1 . The stiffness matrix K stores uncoupled blocks K bd and K c corresponding to each structure, and a coupling term K detailed in section 2.3 reflecting the stiffness of the bearings, denoted B1 and B2, and the fan frame depicted in Fig. 3(b). Finally, f ext stores the prescribed external loads acting on the casing and bladed-disk respectively, and f cn stores the contact forces which couple the two underlying linear systems through the unilateral and frictional contact constraints matrices B.

Craig-Bampton model reduction

In order to gain computational efficiency while capturing the essential dynamic response, the bladed-disk dynamics is reduced through the Craig-Bampton component mode synthesis method [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF]. The DoF x bd are partitioned into internal DoF x i and boundary DoF x b : the latter is kept in the reduced model for contact treatment purposes [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF] and the former is reduced to Á of modal participations q cb using:

x bd D  x b x i à D Ä I 0 ‰ s ‰ c  x b q cb à (3)
where ‰ s stores the static modes. Its size is the number of boundary DoF kept in the reduced model. Matrix ‰ c stores Á constraint modes to control the accuracy of the reduction basis. Both radial and tangential DoF of every blade-tip are kept in the reduced model as well as the two DoF from the center node. Also, Á D 110 constraint modes are incorporated as stated in the next section, leading to 156 DoF for the structural matrices related to the bladed-disk2 of Eq. (1). As for the casing, it is not reduced in this investigation and all 40 nodes are considered for contact treatment.

Bearings and frame coupling methodology

Contrary to the assumptions made in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF], where the bladed-disk and the casing are independent structures solely coupled by the unilateral contact constraints, as displayed in Fig. 3(b), the structural couplings provided by the fan frame and bearings is accounted for in the present study. These two components are incorporated in the model through an array of linear springs, which differs from the approach in [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF] where they are represented by a pair of simple isotropic springs. As detailed in the following, the added complexity in the frame/bearings model introduces time-dependent periodic terms in the stiffness matrix due to the rotation.

As detailed in next section, 1-nodal diameter bladed-disk modes are the only modes potentially exhibiting off-axis shaft motions and participating in the coupling with the casing through the bearings. Hence, it is proposed to explicitly derive the coupling matrix K of Eq. ( 1) in a reduced-order modal space, constraining the bladed-disk distortions to a specific set of modes: the first 1nodal diameter flexural mode .u c ; u s / pictured in Fig. 6(b) and the suspension modes .u c o ; u s o /, both featuring non-negligible shaft motions. For the casing, all nodal diameters actually participate in the coupling and should be incorporated to obtain a general version of matrix K . However, as a first approximation, it is proposed to solely consider the participation of its first 1-nodal diameter mode .v c ; v s / which was shown to be dominant in the interactions [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF]. Accordingly, the physical displacements are projected onto the modal basis through the change of variable x D Vu, where V contains the uncoupled eigenvectors V bd and V c , and u are the corresponding modal coordinates, which in the reduced-order modal basis take the form:

x D Q V Q u (4) with Q u T D v c ; v s ; u c ; u s ; u c o ; u s o and Q V D Ä Q V c 0 0 Q V bd
where Q V contains the associated set of eigenvectors and is of size n 6. Thus, the dynamics of the bladed-disk/casing system is described solely considering modes whose radial distortions contain an important off-axis component at the shaft/bearing interface and may be visualized as the two elastic rings illustrated in Fig. 4, where the blue springs represent the frame/bearings coupling. As adopted in [START_REF] Lesaffre | Contact analysis of a flexible bladed-rotor[END_REF], this choice results in the following displacement fields:

u bd .˛; t / D u c .t / C u c o .t/ cos ˛C u s .t/ C u s o .t/ sin vc .Â; t/ D v c .t / cos  C v s .t/ sin  (5)
where u bd .˛; t / and v c .Â; t/ correspond to the displacements of the bladed-disk and the casing respectively, described by the angles  2 OE0 I 2 and ˛D  t in the static frame of reference. The gap between the two structures g.Â; t/ is expressed as:

g.Â; t/ D ı C v c .t / cos  C v s .t/ sin  u c .t / C u c o .t/ cos. t / u s .t / C u s o .t / sin. t / (6)
where ı is the initial operating clearance introduced earlier. Further, the rotational speed of the bladed-disk is assumed to be constant in counter-clockwise direction and the springs coupling stiffness is denoted Ä. From the assumed displacement forms in Eq. ( 5), Hamilton's principle is adopted to derive the associated equations of motion. The kinetic and strain energies of the two rings are:

T bd .t/ D 1 2 m u P u c .t/ 2 C P u s .t/ 2 C 1 2 m o P u c o .t/ 2 C P u s o .t / 2 T c .t/ D 1 2 m v P v c .t/ 2 C P v s .t/ 2 K bd .t/ D 1 2 k u u c .t/ 2 C u s .t/ 2 C 1 2 k o u c o .t/ 2 C u s o .t/ 2 K c .t/ D 1 2 k v v c .t/ 2 C v s .t/ 2 (7)
where .m u ; m o ; k u ; k o / and .m v ; k v / are the modal masses and stiffnesses of the bladed-disk and the casing, respectively. The strain energy of the coupling stiffnesses Ä arises as:

K Ä .t / D 1 2 Ä Z 2 0 v c .t / cos  C v s .t / sin  u c .t / C u c o .t / cos. t / u s .t / C u s o .t / sin. t / 2 d (8)
Accordingly, the following stiffness is obtained in the modal coordinates:

Q K. ; t/ D Q K C Q K . ; t / D Q V T K Q V C Q K . ; t / D 2 6 6 6 6 6 4
! 2 v C v 0 u cˇ u sˇ u cˇ u s0 ! 2 v C v u sˇ u cˇ u sˇ u cˇ v cˇ v sˇ! 2 u C u 0 0 0 v sˇ v cˇ0 ! 2 u C u 0 0 v cˇ v sˇ0 0 ! 2 o C u 0 v sˇ v cˇ0 0 0 ! 2 o C u 3 7 7 7 7 7 5 (9)
with cˇD cos t and sˇD sin t,

! 2 u D k u =m u , ! 2 o D k o =m o and ! 2 v D k v =m v .
The normalized coupling terms 3 become u D Ä=m u D Ä=m o and v D Ä=m v , and the resulting stiffness matrix-symmetric only under the condition that m v D m u D m o -has time-periodic coefficients. This matrix is the sum of the uncoupled stiffness matrices projected onto the reduced-order modal basis and the frame/bearings coupling matrix, which contains the time-dependent terms. This results in the following reduced equations:

R Q u C Q D P Q u C Q K. ; t / Q u D Q V T .f ext C Bf cn / (10)
In Eq. ( 10), the forcing terms are kept in the physical space for readability: this is where the blade-to-casing clearances and contact forces can be determined precisely. The reduced damping matrix

Q D D Q V T D Q V is diagonal
and the time-dependent stiffness matrix (9) must be computed at each time-step of the time-marching strategy implemented in the following.

Modal analysis of the bladed-disk/casing system 3.1 Shaft motions in cyclically symmetric structures

Perfectly tuned cyclically symmetric structures are those composed of a repetition of N identical sectors. Such structure feature modes that can be sorted by modal families and nodal diameters k-also called spatial harmonics-depending on their frequency and geometrical shape [START_REF] Thomas | Dynamics of rotationally periodic structures[END_REF]. In order to determine which spatial P harmonics are affected by an off-axis motion4 of a node belonging to the axis of rotation of this type of structure, the physical displacements x, limited to the components associated to such node as illustrated in Fig. 5, are translated into the Fourier space O u by means of the change of variable proposed in [START_REF] Bladh | Efficient predictions of the vibratory response of mistuned bladed disks by reduced order modeling[END_REF]. In the local frame of the n-th sector .R n /, this displacement is:

x R n D Â a cos.n˛/ a sin.n˛/ Ã with n D OE1I : : : I N ( 11
)
where a is the displacement amplitude and ˛D 2 =N corresponds to the so-called fundamental interblade phase shift. By manipulating the expressions of the resulting cyclic displacements for each spatial harmonic O u k , it can be shown that: Single harmonics Either k D 0 or k D N=2 (for N even), the cyclic displacements correspond to:

O u 0 D 1 p N N X nD1 x R n D 0 O u N 2 D 1 p N N X nD1 . 1/ n 1 x R n D 0 (12)
which do not participate in the off-axis shaft displacements. Double harmonics Ranging from k D 1 to k D K D bN=2 1c, the associated cyclic displacements are:

k D 2; : : : ; K; O u k c D O u k s D 0 k D 1; O u k c D O u k s D a p 2N cos.˛/ ( 13
)
where the only non-zero contribution corresponds to the first spatial harmonic. It is then clear that a point lying on the rotational axis of a cyclically symmetric structure displays off-axis displacements on 1-nodal diameter modes only. Accordingly, any precessional motion of the shaft is completely described through a combination of 1-nodal diameter modes. Similar conclusions hold for circular plates [START_REF] Kim | Exact solutions of in-plane natural vibration of a circular plate with outer edge restrained elastically[END_REF].

Modal analysis of the uncoupled structures

The modal properties of the bladed-disk and casing are now studied independently. The prediction of modal coincidence in the sense of [START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF] and the corresponding whirling motions involved in the interactions are discussed. For the sake of confidentiality, all frequencies are normalized with respect to the bladed-disk lowest eigenfrequency.

Bladed-disk

The sensitivity of the structure mode shapes and eigenfrequencies to the flexibility of the shaft is now assessed. In agreement with the developments presented above, Figs. 6 display 7 show a comparison of the eigenfrequencies with a free and clamped center node. Clearly, such conditions only affect the 1-nodal diameter modes of each modal family. In conclusion, accounting for the flexibility of the shaft has two main consequences: it generates a shift in 1-nodal diameter eigenfrequencies as these modes can now reveal a whirling component, and it engenders the suspension modes, similar in shape to the k D 1 (1F) mode but featuring a weaker distortion of the blades.

Modal coincidence predictions and whirl motions

Modal coincidence occurs when traveling waves of matching nodal diameters exhibit the same propagation speed [START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF]. The high linear blade-tip velocities make it possible to assume that sliding only occurs during contact phases. Accordingly, the waves stemming from the interactions are always counter-rotating on the rotor and co-rotating on the stator. Therefore, the linear modal coincidence prediction of critical speeds is [START_REF] Batailly | Unilateral contact induced blade/casing vibratory interactions in impellers: Analysis for rigid casings[END_REF]:

! c k D k cr k ! bd k ( 14
)
where ! c k and ! bd k are the casing and bladed-disk eigenfrequencies of matching nodal diameter k, and the corresponding critical speeds cr k are represented by the crossing of modal lines in the Campbell diagram 8.

For the models considered in this investigation, in Fig. 8 several crossings occur in the rotational speed range 2 OE0:7 I 1:3 for the first family of modes, while for the second family, these crossings are spread over a larger frequency range and occur at higher rotational speeds. With regard to interactions between 1-nodal diameter modes that potentially lead to high-amplitude whirling motions, as established in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF], the first flexural (1F) modal family crossing appears at cr 1 ' 1:05, while for the suspension modes it occurs at cr susp ' 2:6, which lies well beyond the operating speed range corresponding to 2 OE0:6 I 2.

According to [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF], depending on the contact location with respect to the shaft displacements, either forward or backward whirling motions are likely to be excited because of the momentum that frictional forces generate at the center of the bladed-disk. In this sense, two scenarios are identified: Standard scenario Contact occurs in the direction of the shaft displacements as illustrated in Fig. 9(a) where the blue arrow indicates the shaft displacement and the blue zone, the preferred contact location. The casing distortion is of smaller magnitude than that of the bladed-disk and the frictional forces (green arrows) create a momentum on the shaft which is in a direction opposite to the rotational speed leading to backward whirling. Inverted scenario In this configuration illustrated in Fig. 10(a), the distortion of the casing is larger than that of the bladeddisk, absorbing the blade-tip clearance in the opposite side of the shaft displacements. This yields forward whirling motions since the momentum generated by the frictional forces is now The particularity of the inverted scenario is that there would be two traveling waves of opposite direction propagating simultaneously on the bladed-disk. Indeed, a co-rotating wave would be generated on the shaft while the frictional contact forces also produce a counter-rotating one on the blade-tips according to the permanent sliding assumption. This suggests that there exists a node of vibration along the radius of the bladed-disk for the change of sign in the traveling waves to occur, hence it is necessary that several 1-nodal diameter modes participate in the response.

For both contact configurations, the kinematic phase angle of the shaft and the shortest blade-to-casing clearance are plotted in Figs. 9(b) and 10(b) over one period of rotation. In the standard scenario, the minimum gap is always in phase with the shaft displacement, while for the inverted scenario, they are also in phase but shifted by an angle D .

Modal analysis of the coupled structures

The modal properties of the coupled structures-through the bearings and frame-described in the modal subspace by Eq. ( 10) are now assessed, thereby neglecting the forcing and damping terms. Hence, the free and conservative system is:

R Q u C Q K. ; t / Q u D 0 (15)
where the time-periodic coefficients 5 of the stiffness matrix defined in Eq. ( 9) can be advantageously transformed into time-independent terms through the change of variable [START_REF] Shapiro | Rotating class of parametric resonance processes in coupled oscillators[END_REF]:

0 B B B B B B @ v c v s u c u s u c o u s o 1 C C C C C C A D 2
q 1 q 2 q 3 q 4 q 5 q 6 1 C C C C C C A (16)
with cˇ0 D cos t 2 and sˇ0 D sin t 2 , corresponding to change of frame at a rotation of half the rotational speed. Therefore, this change of variable implies:

Q u D P.t /q P Q u D P P.t /q C P.t / P q R Q u D R P.t /q C 2 P P.t / P q C P.t / R q (17)
which leads to the following equation of motion with constant coefficients: P q 1 P q 2 P q 3 P q 4 P q

0 B B B B B B @ R q 1 R q 2 R q 3 R q 4 R q 5 R q 6 1 C C C C C C A C 2
5 P q 6 1 C C C C C C A C 2 6 6 6 6 6 6 4 o 0 0 0 0 o 0 0 0 u 0 0 0 0 u 0 0 0 v 0 0 0 0 v 3 7 7 7 7 7 7 5 0 B B B B B B @ q 1 q 2 q 3 q 4 q 5 q 6 1 C C C C C C A D 0 B B B B B B @ 0 0 0 0 0 0 1 C C C C C C A (18) with, o D ! 2 o C 2 4 ; u D ! 2 u C 2 4 ; v D ! 2 v C 2 4
, and D u D v is the normalized coupling stiffness considering that m u D m v D m o . Thus, Eq. ( 18) is written in a compact form as:

R q C G. / P q C K. /q D 0 (19)
corresponding to a set of second-order linear ODE with constant coefficients with:

P 1 P D I 6 G. / D 2P 1 P P K. / D P 1 R P C P 1 Q K. ; t /P (20)
and for which a stability analysis of the equilibrium solution q D 0 becomes straightforward [START_REF] Meirovitch | Computational methods in structural dynamics[END_REF].

Linear stability analysis and mode types

First, Eq. ( 19) is expressed in its state-space form:

P z.t / C Az D 0 (21)
where

z D  q P q à and A D Ä 0 I K. / G. /
The stability of the fixed-point is then addressed through Floquet theory by computing the eigenvalues of the transformation matrix A in Eq. [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF]. For a pair . ; Ä/, the equilibrium position becomes unstable when one of the eigenvalues show a strictly positive real part.

The stability chart in Fig. 11 is computed with respect to Ä and , as the modal properties of the bladed-disk and casing are not considered as variables: unstable responses are represented in red, the color being proportional to the magnitude of the associated eigenvalue real part. Two unstable zones denoted A and B correspond to the coupling of the casing modes with the 1-nodal diameter bladed-disk and suspension modes respectively, where the range of the associated unstable rotational speeds increases with Ä. As detailed in [START_REF] Shapiro | Rotating class of parametric resonance processes in coupled oscillators[END_REF], this instability is a parametric res- onance which depends on the and Ä parameters considered, as well as the eigenvalues of both bladed-disk and casing. It is worth noticing that co-rotating modes do not yield unstable configurations, as both unstable zones emanate from the linear modal coincidence condition ( 14) calculated in the previous section . cr 1 ' 1:05 and cr susp ' 2:6/. The stability diagram 11 suggests that structural changes either on the bladed-disk or on the casing can shift cr 1 out of the operating range in order to avoid instabilities. The stiffness of the frame and bearings between the two structures not only appears to modify the band width of the predicted unstable speeds but it also introduces a shift in the first unstable configuration. The imaginary and real parts of the eigenvalues corresponding to the blue line in Fig. 11 are depicted in Fig. 12. For a given pair .Ä; / leading to unstable motions, three types of modes actually coexist: oscillatory, damped and divergent modes.

Coupled modes and contact initiation

To obtain a visual representation of the bladed-disk/casing system coupled modes, the resulting modal displacements in the timeinvariant modal space must be recast into the original physical space. This is of particular interest for the appropriate characterization of the contact initiation: direction of contact with respect to the shaft displacement and types of waves traveling along the casing and the bladed-disk. The following steps summarize the projection procedure: 1. Complex modal analysis of A in Eq. ( 21) for a pair .Ä; /.

2. Dynamics of a selected mode in time with appropriate initial conditions. 3. Projection of the modal solution into associated Q u coordinates using Eq. ( 16). 4. Projection into the initial physical space x through modeshapes Q V in Eq. ( 4) to obtain the corresponding bladed-disk and casing modal displacements. The mode considered in the projection is the divergent one, arising from the selected .Ä; / parameters that render the equilibrium solution unstable, which are represented by the blue circle in Fig. 11. The associated time-response is projected into physical coordinates through the above procedure. As expected from the stability analysis, this divergent mode is characterized by a backward traveling wave on the bladed-disk and a forward rotating one on the casing, whose deformed shapes are displayed in Fig. 13 and are similar to those of the uncoupled system illustrated in Fig. 9(a).

Indeed, in Fig. 14(c), the comparison between the shaft phase and the smallest operating clearance phase shows how the two structures vibrate in phase and how contact takes place in the same side of the shaft displacement. This contact configuration is in agreement with the kinematic predictions of linear modal coincidence given in Figs. 9(b) for the standard contact scenario. The corresponding divergent whirling orbits are depicted in Figs. 14(a) and 14(b) where the blue arrow indicates the whirling direction. It is shown how a (growing) backward traveling wave in the rotating frame translates into a (growing) forward traveling wave in the static frame because of the bladed-disk rotation. Therefore, instabilities induced by the frame/bearing flexibility might initiate interactions between the casing and bladed-disk, which potentially develop into modal interactions in the sense of [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF][START_REF] Schmiechen | Travelling Wave Speed Coincidence[END_REF].

Contact simulations

Various mechanisms affect the blade-to-casing clearances in service [START_REF] Lattime | High-pressure-turbine clearance control systems: current practices and future directions[END_REF]. These are classified into axisymmetric and asymmetric clearance changes, the former being caused by uniform loads on both casing and rotor (e.g. centrifugal loads) and the latter by non-uniform loads generally acting on the stationary components (e.g. casing ovalization due to thermal gradients, in-take flow aerodynamic loads). In this context, two contact configurations, both leading to asymmetric clearance changes and potentially producing harmful interactions are now investigated: (1) initial casing distortion, where an external static load is applied on the casing and (2) mass imbalance on the bladed-disk.

Solution algorithm

The numerical procedure implemented in this investigation can be described by the following steps: 1. A perturbation is applied to the system to absorb the initial blade-tip clearance ı and initiate contact, either:

The casing is distorted along a 1-nodal diameter during a short time interval; the two structures are then left free to interact.

A mass imbalance force is applied to the bladed-disk.

2. Equations of motion are solved through an explicit time-marching scheme [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF] and can be treated in two bases:

In the physical space through Eq. ( 1) as in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF].

In the truncated modal space through Eq. ( 10), where the proposed formulation of frame/bearings coupling can be incorporated. 3. If the frame/bearings coupling is included in the formulation, the stiffness matrix in Eq. ( 9) is calculated at each time-step. 4. When a penetration is detected, the contact forces f cn are computed using a Lagrange multipliers method in the physical space [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]. 5. Friction is accounted for through a Coulomb law and sliding is assumed. Based on previous work [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF], the solution algorithm 1 is used to compute the time response of the system.

Space-time discretization convergence

In order to ensure the accuracy of the results, a convergence study was carried out with regards to the spatial and temporal discretization. In this sense, two parameters are of interest: the size of the reduced-order model Á and the time-step t in the time integration procedure.

The nonlinearity of the contact force and the explicit nature of the time-stepping scheme make the algorithm strongly sensitive to the time-step size. Convergence is reached in displacements and contact forces for t D 10 6 s which is preserved for all simulations performed in this investigation.

Similarly, convergence is verified for increasing Á. As in [START_REF] Batailly | Assessment of reduced models for the detection of modal interaction through rotor stator contacts[END_REF], where the notion of motion convergence is introduced, asymptotic convergence is not perfectly achieved in displacements but is reached in the types of motion detected even for a small reduction basis, since the lower frequency modes have a major participation

Results and discussion

The nonlinear dynamics of the bladed-disk/casing system is first investigated in the full space (Eq. ( 2)) with respect to the abovementioned perturbation mechanisms and where the frame/bearings coupling is omitted as in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF]. The system response is then considered in the reduced-order modal basis through Eq. [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF], where the effects of the kinematic constraints imposed by the modal projection are discussed and the role of the frame/bearings coupling is addressed.

Scenario A: contact initiation through casing distortion

The first scenario involves a casing statically distorted into a knodal diameter shape in order to absorb the initial blade-tip clearance [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF]. A casing ovalization .k D 2/ generates a symmetric response on the rotor with two diametrically opposed blades ex-hibiting an identical response: this results into a zero displacement of the center node and is not further discussed. Instead, it is assumed that the casing distortion is along a 1-nodal diameter shape, which may be caused by aerodynamic loads at take-off [START_REF] Lattime | High-pressure-turbine clearance control systems: current practices and future directions[END_REF] and which seems to favor high-amplitude whirling motions as reported in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF] where three motion types are predicted: The simulations are carried out considering the parameters listed in Tab. 2, where the structural coupling Ä is set to zero. The resulting frequency response of the bladed-disk is depicted in Fig. 15, evaluating the center node displacements .F x o /. The FFT is computed once the transient response is damped out. This figure features a first group of critical speeds around the linearly predicted one cr 1 ' 1:05 for which the highest amplitudes are reached. In general, it is shown that the dynamics is mainly driven by the first family of modes regardless of the rotational speed, with a minor contribution of the suspension modes and a negligible participation of all other modes. Other critical interactions are also detected-denoted cr nl in Fig. 15-that cannot be predicted by the linear modal coincidence criterion.

Divergent
The linear modal interaction at cr 1 is further analyzed in Fig. 16 where the associated whirling orbits of the shaft are depicted. The interaction is characterized by a backward whirl in the rotating frame and the whirling direction is unchanged in the static frame. However, since the first eigenfrequency of the casing ! c

1 is very low and the resulting critical speed cr 1 is very close to ! bd 1 , the precessional motion of the shaft is almost inexistent in the static frame as shown in Fig. 16(b). The spectrogram of the response is displayed in Fig. 17(a) along with the modal strain energy associated to the first family of modes, the suspension modes and the global response of the bladed-disk in Fig. 17(b). This energy is calculated as:

SE bd D 1 2 x bd T Kx bd D 1 2 u bd T ƒu bd with ƒ D V T KV (22)
where u bd contains the modal contributions of specific modes 6 and ƒ the eigenvalues associated to the eigenmodes V of the system. Spectrogram 17(a) shows how the suspension modes have a non-negligible contribution during the transient part of the response t 2 OE0 s I 0:5 s while the remaining of the response lies on the 1F modal family. This is consistent with the energy distributions in Fig. 17(b) where it is clear that the 1-nodal diameter mode is responsible for the divergence (as expected), with a very rapid growth in amplitude after t D 1:8 s. It is also shown that the k D 0 mode is present in the divergent part of the response along with the suspension modes, while all other nodal diameters do not seem to participate. However, during the transient part of the response higher order modes also participate in the response, as it may be seen that the 1F modal family and the suspension modes cannot account for the total strain energy.

Further post-processing of the other critical interactions mentioned above result in a very similar behavior in terms of strain energies, vibration amplitudes and frequency content, not shown here for the sake of brevity. However, as depicted in Figs. 18(a) and 18(b) for cr nl D 1:39, a backward whirling orbit is obtained in the rotating frame similar to the case of cr 1 , which is translated into a forward traveling one in the static frame due to the higher rotational speed. In this particular case, the 1-nodal diameter of 6 All modes are considered when calculating the global strain energy the 1F modal family of the bladed-disk interacts with the 5thsuper-harmonic in time of the mode of casing .! c D 5! c 1 /. This type of interaction cannot be predicted by the linear modal coincidence criterion ( 14) but exhibits a similar order of magnitude in terms of displacements than those detected for cr 1 , hence being just as critical from a structural integrity perspective.

In terms of contact location versus shaft displacement, it is shown in the phase comparison depicted in Fig. 18(c) that the system responds as in the standard contact scenario described in the previous section. Although in this case, for the mode responsible of the divergence whose modeshape is illustrated in Fig. 6(b), the most elongated blade is not perfectly aligned with the shaft maximum radial displacement, thus introducing the phase shift between the two quantities visible in Fig. 18(c). A decrease in the whirling frequency is also visible in this figure: there are three oscillations for t 2 OE2 s I 2:2 s while there are only two for t 2 OE2:8 s I 3 s. This behavior can be explained by the frequency increase of the bladed-disk response in the relative frame during the divergent phase of the simulation, visible in spectrogram 17(a) for D 1:05, which translates into a frequency decrease in the fixed frame due to the wave propagation direction. A similar behavior was observed for all the other divergent interactions detected. Further, analogous conclusions than in [START_REF] Salvat | Two-dimensional modeling of shaft precessional motions induced by unilateral blade/casing contacts in aircraft engines[END_REF] can be drawn for other types of regimes, as sustained behaviors are observed near critical speeds and damped motions are attained elsewhere in the rotational speed range. It should also be noted that the shaft dynamics appeared to be crucial in the emergence of divergent regimes, as a clamped center node configuration led only to damped motions with the same simulation parameters listed in Tab. 2.

Scenario B: contact initiation through mass imbalance

The mass imbalance is reflected via an external forcing term with a constant component in the rotating frame and proportional to the square of the rotational speed. Its magnitude is selected so that the initial blade-tip clearance is absorbed with the lowest value considered in Tab. 2. The main difference with scenario A is that the zero solution of Eq. ( 1) no longer exists and therefore, the former damped regimes cannot be reached. This type of contact initiation mechanism is only meaningful when extending the scope of the interactions to whirl motions and hence could not be studied with a rigid shaft assumption considered in [START_REF] Legrand | Two-dimensional modeling of an aircraft engine structural bladed disk-casing modal interaction[END_REF]. The same post-processing strategy as for scenario A is considered. The frequency response of the shaft is depicted in Fig. 19 over the entire rotational speed range without accounting for the frame/bearings coupling but with the input parameters listed in Tab. 2. In this figure, the same colormap as in Fig. 15 is used in order to highlight the change in vibration amplitudes, which are generally smaller than those obtained when a modal interaction arises. The shaft is shown to respond mainly on the suspension modes and the 1F modal family, and some interaction regimes are detected for 2 OE1:42 I 1:75 where the amplitudes become of the same order of magnitude as in the first configuration.

For the critical speed cr 1 , the shaft displacements in the rotating frame are depicted in Fig. 20(a) along with whirling orbit in the static frame in Fig. 20(b). The shaft whirls in an circular X o / with a sampling frequency equal to the rotational speed D 1:05, where the limit cycle oscillation of the steady-state response evolve in a highly localized zone of the .X o ; P X o / plane, thus suggesting that the shaft exhibits an almost perfectly periodic motion. In fact, the strain energy distribution depicted in Fig. 20(e) shows how the mass imbalance load mainly excites the suspension modes while the remaining of the strain energy is localized on the 1-nodal diameter mode, all other modes having a negligible contribution to the overall dynamics. It is expected for the mass imbalance to mainly excite the suspension modes in a linear framework, yet in the nonlinear framework introduced by the unilateral contact constraints, all modes of both structures tend to be excited by the contact loads. For the current models, at the predicted modal coincidence speed, the arising rotor/stator interactions do not seem to develop into divergent regimes and are well contained by the surrounding casing (vibration levels are in the order of 10 µm for the casing). It appears that the energy levels stabilize after t ' 1:1 s as steady-state is reached and are about ten times lower in comparison to Fig. 17(b), where for the same rotational speed, the transient excitation develops into a modal interaction resulting into a continuous growth of vibration amplitudes.

The blade-to-casing distances are displayed in Figs. 21(a) and 21(b) where it is shown that: several blades remain in permanent contact; two blades exhibit intermittent contacts and the remaining ones vibrate freely, solely excited by the disk coupling with the neighboring blades. This behavior is consistent with the findings reported in [START_REF] Lesaffre | Contact analysis of a flexible bladed-rotor[END_REF][START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF], where the resulting contact loads are almost constant on the bladed-disk and act as rotating loads on the casing. However, in [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF], these operating conditions appeared to induce a growth of vibration amplitudes on both structures at the modal coincidence speed, combining the whirling motion of the shaft with a 2-nodal diameter wave on the casing and the blades. This major difference with respect to the behavior observed in the current investigation might be explained by the rigid disk assumption and the clamped boundary conditions imposed at the blades foot in [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF], and as suggested in [START_REF] Batailly | Assessment of reduced models for the detection of modal interaction through rotor stator contacts[END_REF], by the strong kinematic constraints induced on the models in [START_REF] Parent | Whole engine interaction in a bladed rotor-to-stator contact[END_REF] which tend to artificially favor the detection of modal interactions. Additionally, Nevertheless, some high-amplitude interaction regimes are indeed detected for 2 OE1:42 I 1:75. It is proposed in the following to characterize the response in speed range and whether or divergent regimes are In the distribution for D 1:6 shown in 22, it clearly visible how the energy levels are much higher than those depicted in Fig. 20(e) for D 1:05 and are in the same order of magnitude than when a modal interaction arises in Fig. 17(b). The main dif-ference between the two types of interactions is that the vibrations are contained by the casing when a mass imbalance is the source of excitation, while the vibrations grow exponentially when modal occurs. Indeed, as displayed in Fig. 22, a rapid growth in the strain energy is detected after t D 0:5 s, reaching a peak in amplitude of about 800 J over the time interval t 2 OE1:8 I 2:2 s and decaying the end of the simulation 7 . As for the modes participating this interaction, it is shown that: besides the contribution of the suspension modes which are dominant during the first 0:2 s of the response and exhibit rather constant energy levels throughout the rest of the simulation, of the k D 0 and k D 1 modes is noted during the transient growth in energy levels as well as a non-negligible participation of all other harmonics throughout the entire response.

Sensitivity to the frame/bearings coupling

Lastly, the nonlinear dynamics of the bladed-disk/casing reducedorder coupled model are investigated with respect to the two loading scenarios discussed above for the full uncoupled system, the main objective being to characterize the influence of the frame/bearings stiffness Ä on the interactions detected.

Casing distortion excitation

The linear modal coincidence captured in the full space at cr 1 depicted in Figs. 16 and20, is also detected in the reduced space as shown in Fig. 23 for Ä D 0 N=m, since the modes kept in the truncation are those having a major participation in the system response. It may be seen in this figure, in agreement with prior observations, how the 1-nodal diameter mode is the one responsible for the divergence. It is followed by a very rapid growth in amplitude of the suspension modes strain energy, which generate a far more aggressive response than with the full system. Indeed, as explained in [START_REF] Batailly | Assessment of reduced models for the detection of modal interaction through rotor stator contacts[END_REF], the kinematic restrictions imposed by the modal projection artificially favor the initiation of high-amplitude interactions.

The sensitivity of the bladed-disk response to the frame/bearings coupling stiffness Ä is analyzed in Fig. 24, where the maximal displacement of the center node is displayed for rotational speeds near cr 1 and three different values of Ä D OE0 I 500 I 1000 N=m. It is shown in this figure how the structural coupling introduces a stiffening of the system and produces a shift in the modal interaction speed, as would be expected from the stability analysis carried out in section 3.3. The frequency content and energy distributions, as well as the contact direction with respect to the shaft displacements, are not affected by the coupling and again the 1-nodal diameter of the 1F family generates the divergences.

Mass imbalance excitation

The coupling matrix seems to have a negligible influence on the system response. As depicted in Fig. 25(a), very similar shaft displacements are obtained at cr 1 for all values of Ä D OE0 I 500 I 1000 N=m tested. Similar than for the full system, these displacements correspond to a forward whirling circular orbit in the static frame of about 1:2 mm radius with oscillations in the order of 0:1 mm. As for the modal deformations illustrated in Fig. 25(c), these are also in agreement with prior observations in terms of energy levels and localization, since the suspension modes are largely dominant throughout the entire response. Also, regarding the interactions encountered in the speed range 2 OE1:42 I 1:75 for scenario B, these could not be observed in the truncated modal space, since they were characterized by a transient growth of vibration amplitudes and a non-negligible participation of all the 1F harmonics, in particular the k D 0 and k D 1 as was shown in Fig. 22. Instead, divergent regimes were obtained for > 1:52 where, as illustrated in Fig. 26, the strain energy distribution and associated magnitudes become similar to those of modal interactions (c.f. figure This indicates that the truncation modal basis chosen in Eq. ( 4) is not sufficiently rich These observations suggest that the frame/bearings coupling between casing and bladed-disk plays a major role in the modal interaction regimes detected in scenario A and in the subsequent high-amplitude precessional motions of the shaft, and has a limited influence in the system behavior when a mass imbalance is the source of excitation. The modal truncation basis proved to be sufficient to capture the dominant dynamics of the system in scenario A. However, even this basis is appropriate in scenario B for small rotational speeds, for > 1:52 non-physical divergences were reached regardless of the coupling stiffness Ä considered due to the strong kinematic constraints imposed by the modal projection.

Conclusions and perspectives

This study focuses on the occurrence of shaft precessional motions induced by unilateral blade/casing contacts in a two-spool commercial engine fan stage. A 2D in-plane FE-model representative of a bladed-disk/casing system is built, where the flexibility of the shaft is represented by a pair of linear springs linked to the center of the disk which drive the so-called suspension modes. Additionally, it is proposed to model the structural link between bladed-disk and casing provided by the fan frame and the bearings through an array of linear springs. The associated methodology, based on the derivation of Hamilton's principle for the calculation of the corresponding coupling stiffness matrix, involves a truncated set of modes utilized for the projection of the equations of motion. This coupling stiffness matrix contains time-periodic coefficients which can be eliminated through a change of frame in order to carry out a linear modal analysis, from which the stability of the zero solution revealed two unstable zones where damped, oscillatory and divergent modes coexist.

By means of a time-marching strategy, the unilateral contact interactions between these flexible structures are initiated via two distinct mechanisms: (1) a prescribed casing distortion and (2) a mass imbalance on the bladed-disk. The shaft dynamics proved to have a key role in producing potentially harmful regimes, in particular in the former scenario which can lead to divergent modal interactions. For most of the interactions detected, the shaft precessional motion takes the form of a backward traveling wave in the rotating frame which is translated into a forward traveling one in the static frame. Through a projection from physical to cyclic coordinates and the calculation of the associated modal strain energy, it is shown how the 1-nodal diameter mode of the first modal family of the bladed-disk is dominant when a modal interaction occurs, while the suspension modes have a major role when a mass imbalance is the source of excitation.

In the truncated modal space, the nonlinear dynamics of the reduced-order coupled system are studied considering both loading scenarios. The modal projection proved to be suitable for the detection of modal interactions initiated by the casing distortion, where the frame/bearings stiffness appeared to introduce a shift in the detected critical speeds in agreement with the linear stability predictions. However, when a mass imbalance is the source of excitation, the obtained results are shown to be valid up to a certain speed . > 1:52/ beyond which only divergent regimes are observed, as the strong kinematic constraints induced by the modal projection do not allow to capture the behavior of the full system where an important contribution of all harmonics is noted.

These qualitative observations set a basis for two different research tracks which are currently being pursued: firstly, one regards an extension of the coupling matrix to all nodal diameter modes in order to avoid the modal truncation and study the dynamics of the fully coupled structures. This will enable to better comprehend the influence of the coupling stiffness in the detected interactions and a generalization of the stability analysis here presented. Geometric nonlinearity effects are also to be included in the formulation. The second track concerns the implementation of 3D industrial bladed-disk/casing models, which exhibit more realistic and complex veering behaviors, that will allow to account for out-of-plane motions and for the whirling conical component induced by the shaft bending. This type of model will also include centrifugal and gyroscopic effects, which are neglected in this two-dimensional study, and the role of the abradable coating wear on the high-amplitude shaft precessional motions arising from the contact events will be investigated.

Figure 1 .

 1 Figure 1. Turbofan engine components: 1. fan; 2. Low Pressure (LP) compressor; 3. High Pressure (HP) compressor; 4. combustion chamber; 5. HP and LP turbines

Figure 3 .

 3 Figure 3. Fan stage model

Figure 4 .

 4 Figure 4. Equivalent reduced-order model: deformed [ ] and undeformed [ ] bladed-disk, deformed [ ] and undeformed [ ] casing and frame/bearings [blue springs]

Figure 5 .

 5 Figure 5. Example of cyclically symmetric structures with 6 sectors

Figure 6 .

 6 Figure 6. First family of free vibration modes: center node displaced [ ] and rest [ ] position

  First two families of modes

Figure 7 .Figure 8 .

 78 Figure 7. Veering diagram with clamped [ ] and free [ ] center node

Figure 9 .

 9 Figure 9. Regular contact configuration for 1-nodal diameter modes

Figure 10 .

 10 Figure 10. Inverted contact configuration for 1-nodal diameter modes

Figure 11 .

 11 Figure 11. Stability analysis: [ ] unstable configurations and [ ] example point

Figure 12 .

 12 Figure 12. Evolution of eigenvalues for Ä D 500 N=m

! bd 1 ! c 1 timeFigure 13 .Figure 14 .

 111314 Figure 13. Deformed 1-nodal diameter mode shape in time

Figure 15 . 1 Figure 16 .

 15116 Figure 15. Frequency response of center node for casing distortion excitation

Figure 17 .Figure 18 .

 1718 Figure 17. Results for initial casing distortion at D 1:05 and Ä D 0 N=m

Figure 19 .

 19 Figure 19. Frequency response of center node for mass imbalance excitation

Figure 20 .

 20 Figure 20. Results for mass imbalance load at D 1:05 and Ä D 0 N=mco-rotating orbit of about 1:1 mm radius with oscillations in the order of 0:1 mm. The associated in Fig.20(c), defined by the lateral shaft displacements X o and velocities P X o , clearly shows that once the transient of the response damps out .t > 1:1 s/ the solution reaches a stable limit cycle. This is further highlighted in the Poincaré map 20(d), defined in the state-space .X o ; P X o / with a sampling frequency equal to the rotational speed D 1:05, where the limit cycle oscillation of the steady-state response evolve in a highly localized zone of the .X o ; P X o / plane, thus suggesting that the shaft exhibits an almost perfectly periodic motion. In fact, the strain energy distribution depicted in Fig.20(e)shows how the mass imbalance load mainly excites the suspension modes while the remaining of the strain energy is localized on the 1-nodal diameter mode, all other modes having a negligible contribution to the overall dynamics. It is expected for the mass imbalance to mainly excite the suspension modes in a linear framework, yet in the nonlinear framework introduced by the unilateral contact constraints, all modes of both structures tend to be excited by the contact loads. For the current models, at the predicted modal coincidence speed, the arising rotor/stator interactions do not seem to develop into divergent regimes and are well contained by the surrounding casing (vibration levels are in the order of

  Blade-to-casing distances: permanent contact [ ], intermittent contact [ ] and free [ ]

Figure 21 .

 21 Figure 21. Contact characteristics: D 1:05 and Ä D 0 N=min the phase comparison depicted in Fig.21(c), it is shown how unilateral contact takes place in the same direction as the shaft displacement once a steady-state is reached. These observations are in agreement with the kinematic predictions given in the previous section in Fig.9(b) and indicate that the imbalance is prone to initiate contact in the expected standard direction.Nevertheless, some high-amplitude interaction regimes are indeed detected for 2 OE1:42 I 1:75. It is proposed in the following to characterize the response in speed range and whether or divergent regimes are In the distribution for D 1:6 shown in 22, it clearly visible how the energy levels are much higher than those depicted in Fig.20(e) for D 1:05 and are in the same order of magnitude than when a modal interaction arises in Fig.17(b). The main dif-

Figure 22 .

 22 Figure 22. Strain energy at D 1:6 and Ä D 0 N=m: total [ ], k D 0 [ ], k D 1 [ ], other k [ ] and susp. [ ]

Figure 23 .Figure 24 .

 2324 Figure 23. Results for casing distortion excitation at D 1:05 and Ä D 0 N=m

Figure 25 .Figure 26 .

 2526 Figure 25. Results for mass imbalance excitation at D 1:05

Table 1 .

 1 Geometry parameters

	Parameter Value bladed-disk radius R bd D 0:5 m disk radius R disk D 0:05 m clearance ı D 1 mm beams cross-section 100 10 mm 2

Table 2 .

 2 motion Vibration amplitudes increase in time, thus involving whirling orbits of growing amplitude: this is the most critical scenario. Sustained motion Intermittent contacts lead to a shaft whirling orbit of almost constant amplitude with a radius of about the initial gap ı. This motion type is generally witnessed for rotational speeds close to critical regimes and the associated vibration amplitudes are sufficiently large to generate highcycle fatigue issues. Damped motion It is characterized by decreasing vibratory amplitudes in time subsequently producing a loss of contact. The interaction is limited to the transient part of the response. Simulation parameters -scenario A

	Parameter Value
	simulation time T tot D 3 s coupling Ä D 0 N=m gap ı D 1 mm friction coefficient D 0:1 rotational speeds 2 OE0:95 I 2 casing distortion k D 1 modal force F c D 30 kN

Modal viscous damping is such that D becomes diagonal in the modal space with non-vanishing terms

! i , where ! i is the undamped natural frequency of the i -th mode and is the damping coefficient.[START_REF] Félix Patrón | New methods of optimization of the flight profiles for performance database-modeled aircraft[END_REF] Without any explicit notice, the structural matrices of the bladed-disk are the reduced version of the finite element ones in the paper.

This section explores the modal properties of the introduced bladeddisk/casing system. First, it is shown that the off-axis displacement components of the shaft are solely described by 1-nodal diameter modes in a cyclic-symmetry framework. Second, a linear modal analysis of the uncoupled structures is carried out and modal coincidence critical speeds are predicted. Last, the dynamics of the coupled structures is analyzed from a modal perspective in order to understand how the frame/bearings coupling affects the overall dynamics.[START_REF] Sobieszczanski-Sobieski | Multidisciplinary aerospace design optimization: survey of recent developments[END_REF] In the projection onto the modal space, the eigenvectors are generally normalized with respect to the mass matrix, hence the corresponding modal masses are unitary and the coupling terms would simply become u D v D Ä.

Only displacements perpendicular to the axis of rotation are considered.

Eq. (15) is a type of Mathieu equation.

These simulations were conducted over larger simulation times and vibration levels stabilize at the same order of magnitude as for D 1:05. For this source of excitation, no divergent regimes were detected.

Acknowledgements

Thanks go to SNECMA for its technical and financial support. This work takes place in the framework of the MAIA mechanical research and technology program sponsored by CNRS, ONERA and SAFRAN Group.