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

The present work targets shaft whirling motions induced by direct blade/casing unilateral contact occurrences in aircraft engine bladed-disk assemblies. These contact events are favored by increasingly reduced blade-tip clearances and potentially lead to harmful interactions that may threaten the engine structural integrity. A simplified 2D in-plane finite element model representative of the engine fan stage is built, accounting for the flexibility of the shaft through two linear springs attached to the disk center node and the structural coupling provided by the fan frame and the bearings, modeled by an array of linear springs. A linear stability analysis of the reduced-order coupled system reveals two unstable zones in a selected rotational speed range, emanating from the linearly predicted modal coincidence speeds. Through a time-marching strategy, two asymmetric contact initiation mechanisms are investigated: (1) a prescribed casing distortion and (2) a mass imbalance on the bladed-disk. It is shown how the 1-nodal diameter mode of the first modal family of the bladed-disk is dominant when a modal interaction arises from the transient casing distortion and leads to divergent regimes. The presence of the frame/bearings coupling induces a shift in the critical speeds detected, generally characterized by a backward traveling wave in the rotating frame and a forward traveling one in the fixed frame. Further, when a mass imbalance is the excitation source, the suspension modes appear to have a major role and a stable limit cycle is reached regardless of the coupling stiffness with much lower energy levels than in divergent regimes.


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) [1], planning cost-efficient flight trajectories [2] as well as the development of multidisciplinary optimization algorithms [3], 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 [4] 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 [5]. As detailed in [6], 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 [7,8].
rotor bearing (a) Rotor/bearing [9] bladed-disk casing (b) Bladed-disk/casing [10] 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 [11] or hydrodynamic [12] bearings, liquid or gas annular seals [13], 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 [8,14]. Generally speaking, three distinct phenomena arising from rotor/stator interactions are discussed in the literature: (1) modal coincidence [15], 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 [16], (2) rubbing, refers to direct contact either between shaft and journal bearings [7] or blade-tips and surrounding casing [17] depending on the system of interest and (3) whirl [9], 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 [7]. 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 [18] 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 [19], 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 [20] 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. [21] 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 [22].
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: (1) an initially prescribed casing distortion and (2) a mass imbalance on the bladed-disk. Based on previous work [23], the interactions between these two flexible structures are studied by means of a time-marching strategy similar to [10]. 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.

Bladed-disk and casing models
The fan stage model in this investigation extends the one introduced in [10] 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 [10], 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 ı.
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 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: that is, in a compact form: 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 [24].
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 [23] and the former is reduced to Á of modal participations q cb using: 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-disk 2 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 [23], 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 [21] 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 [23]. 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: 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 [20], this choice results in the following displacement fields: 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: 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: 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: Accordingly, the following stiffness is obtained in the modal coordinates: with cˇD cos t and sˇD sin t, 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: 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 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.

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 [25]. In order to determine which spatial  harmonics are affected by an off-axis motion 4 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 [26]. In the local frame of the n-th sector .R n /, this displacement is: 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: 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; : : : 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 [27].

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 [15] 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. 4 Only displacements perpendicular to the axis of rotation are considered.

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 a center node displacement for k D 1 only. Additionally, the two DoF of the center node drive the suspension modes through the linear springs connected to it. These modes arise between the first two modal families (green dashed line in Fig. 7(b)) with a frequency f o ' 2:5 and have a shape similar to Fig. 6(b), thus exhibiting a major whirling component. The veering diagrams in Fig. 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 [15]. 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 [28]: 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 [23], 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 [21], 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: 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 [29]: 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: which leads to the following equation of motion with constant coefficients: , 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: corresponding to a set of second-order linear ODE with constant coefficients with: and for which a stability analysis of the equilibrium solution q D 0 becomes straightforward [30].

Linear stability analysis and mode types
First, Eq. (19) is expressed in its state-space form: where The stability of the fixed-point is then addressed through Floquet theory by computing the eigenvalues of the transformation matrix A in Eq. (21). 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 [29], 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  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 [10,15].

Contact simulations
Various mechanisms affect the blade-to-casing clearances in service [6]. 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 [10] and can be treated in two bases: In the physical space through Eq. (1) as in [23].
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 [31]. 5. Friction is accounted for through a Coulomb law and sliding is assumed. Based on previous work [10], the solution algorithm 1 is used to compute the time response of the system.

Input:
geometrical and mechanical properties of casing and bladed-disk reduction basis: Á coupling stiffness: Ä simulation parameters: speed range, time-step, simulation time T tot , perturbation for contact initiation, either: -casing distortion: k, amplitude, loading duration mass imbalance: amplitude initial conditions

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 [32], 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 [23]. The system response is then considered in the reduced-order modal basis through Eq. (10), 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 [23]. 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 [6] and which seems to favor high-amplitude whirling motions as reported in [23] where three motion types are predicted: Divergent 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.  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.
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  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: 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 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 [23] 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 [10].  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 the whirling orbit in the static frame in Fig. 20(b). The shaft whirls in an circular  Fig. 20(c), defined by the lateral shaft displacements X o and velocities P X o , clearly shows that once the transient part 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 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 [20,21], where the resulting contact loads are almost constant on the bladed-disk and act as rotating loads on the casing. However, in [21], 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 [21], and as suggested in [32], by the strong kinematic constraints induced on the models in [21] which tend to artificially favor the detection of modal interactions. Additionally,  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 further characterize the system response in this speed range and determine whether or not divergent regimes are observed. In the energy distribution for D 1:6 shown in Fig. 22, it is 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 coincidence 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 to a minimum by 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, a major influence 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 and 20, 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 [32], 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 23(c)). 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 if 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.