Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques 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
Losing a blade from a turbofan rotor of a jet engine during its normal operation either due to bird-strike or due to fatigue damage is a real possibility, and as such it is a major safety concern in commercial civil aviation [1][2][3][4][5][6]. This possibility has led to the certification requirements of demonstrating the engine capability by running a successful Fan Blade-Out (FBO) test and being able to shut-it off successfully in 15 s. Recently, due to certain field events this requirement has been made more stringent by adding the clause on analyzing the engine blade-out loads during the windmilling of a damaged rotor [7] as a consequence of the fan blade-out event. The civil aviation regulatory agencies such as US Federal Aviation Administration has additional regulation (CFR 33.74) requiring the engine manufacturer to show that continued rotation of the rotors after an FBO event for the remaining diversion mission (up to 3 or more hours) will also result in a safe condition. Other international regulatory agencies such as European Aviation Safety Agency (EASA) also have similar requirements. FAA [8] advisory circular, AC33.74-92-1A, provides guidance on how to demonstrate compliance with the continued rotation regulation. This demonstration is normally done by analysis.
Loss of a fan blade from a turbofan rotor introduces many different types of nonlinearities in the dynamical system. First, the release blade applies a large impact load on the containment system, which is generally a ring type of structure. During this phase of the impact event, depending upon the strength of the containment structure the release blade may or may not be able to punch a hole in the inner shell of the fan casing [9,10]. Most of the test data related to the full fan

Mathematical formulations
In the current paper, we consider dynamic response of an overhung rotor supported by a set of N c number of bearings at multiple axial locations and deformable blades rubbing against a case at the other end (see Fig. 1). The rotor shaft may be under a periodic sinusoidal axial force of magnitude P t superimposed on a constant static axial load P 0 and torque T 0 on the system. In rotating machinery, such pulsating axial forces are usually caused by misalignment, aerodynamic forces, and inlet distortion etc. The rotor system considered in this paper is composed of a flexible hollow shaft, a rigid disk, and a full set of elastically deformable blades, a sector of which may be rubbing against the inside surface of an outer cylindrical containment structure. In this study, we have assumed that the containment structure behaves like a rigid wall and the spinning shaft-disk with bladed rotor as an impactor hitting its inside concave surface at a shallow oblique angle. The dynamics of rotating blade-tips hitting a rigid cylindrical wall was considered in full detail in author's previous investigations [45,46]. The inner surface of the fan case has a lining of a relatively soft material like honeycomb, foam or, an abradable coating against which the blade tip is allowed to rub with stiffness in the radial direction represented as K Rub . The tip-rub forces can deform the blade in all the three directions of the fan-rotor viz., tangential, axial and radial. The radial force at the blade-tip F Rub calculated in this fashion would also result into a tangential load of magnitude mF Rub acting in the direction opposite to the direction of shaft rotation. The effective coefficient of friction m at the blade-tip/ abradable contact interface surface is usually established in an empirical way by measuring the instantaneous speed-drop of the rotor-shaft during the rub-event.
It should be recognized that in determining the dynamic characteristics of an overhung rotor-shaft, the gyroscopic moments play a significant role. This effect is further magnified in an asymmetric rotor with a large-size fan blade missing from one of the disk-sectors caused by a blade-loss. Appendix A describes the energy method approach to derive the mathematical relationship to determine the gyroscopic effects in a bladed nonsymmetric fan-rotor. As it can be seen that this effect not only makes the mass matrix time-dependent varying at twice the frequency of running speed rather gyroscopic matrix also becomes a sinusoidal function of time.
The Coulomb friction between the blade tip and the outer surface introduces skew-symmetric terms in the stiffness matrix for the rotor shaft which, depending upon the amount of effective damping present in the dynamical system, can initiate instability in the bladed-rotor shaft. The basic model has been described in the author's first paper of 2004 on the related topic [31]. In the current paper, the shaft is assumed to be supported at multiple axial locations z i along the length of the shaft ', which may be either ball-bearing or roller-bearing, with their respective stiffness S ðiÞ b and damping D ðiÞ b parameters. As a general case, these support bearing may not necessarily be symmetrical in the X-Y plane and as such for any typical ith bearing located at Z¼ z i their effective stiffness and damping in the matrix form are written as, where, matrix terms such as S bðiÞ xx , S bðiÞ yy , S bðiÞ xy , and S bðiÞ yx have the dimensions of linear stiffness as force per unit length, and B bðiÞ xx , and B bðiÞ yy have the dimensions of bending stiffness as moment/unit radian. For, roller bearings with no rotational resistance the terms B bðiÞ xx and B bðiÞ yy will be zero, but for ball bearings with rotational load carrying capability the terms B bðiÞ xx and B bðiÞ yy are nonzero. VanDuyn [47] discusses the current design approach and its stiffness and damping characteristics for such a bearing support arrangement used in a typical turbofan rotor for specific applications to an aeroengine. 2.1. Equations of motion of a rotor-shaft with deformable pre-twisted blades In the current derivation, we assume that the turbofan rotor is in the form of an overhung cantilever beam with N number of flexible blades mounted on a rigid disk at Z ¼ '. Gennaretti and Muro [48] have used a reduced order formulation for a multibladed rotor system which takes into account of aeroelastic deformation in a rotating frame of reference. Thus, with the usual shaft and bearing notations, we have following two equations to represent the shaft motions in the global X (horizontal) and Y (vertical) directions respectively: external forces coming from the clamped roots of individual blades in the lateral and longitudinal directions for each blades are, À[F Z ] s ¼ 0 and À [F z ] s ¼ 0 , respectively. Here again, the blade-count parameter j varies from 1 to N in a fullybladed symmetric rotor, and the limit changes to (NÀ 1) for a bladed-disk with a missing blade after a fan blade-out event for an asymmetric rotor. In this formulation, the blade airfoils on the fan-disk are modeled as pre-twisted thin shells [44], which are cantilevered with the clamped-end attached to the top surface of the disk at radius 'r'. The actual turbofan blade will have a total pre-twist of (b L À b 0 ) radians with a negative sign. It should be noted that the sample blade shown in Fig. 1 has a positive b, which is consistent with the right-hand vector rotation rule. These blades can bend both along the span and chord as well as twist due to transverse deformation Z(s,c,t) in the thickness direction and can also have membrane stretching due to longitudinal deformation z(s,t) in the span direction of the airfoil. Thus, for an axis-symmetric rotor, we will have N sets of following two equations representing each individual fan blade deformation in their own local lateral and longitudinal directions: and, In the above Eqs. (6) and (7), the force terms [F Z ] s ¼ 0 and [F z ] s ¼ 0 represent pseudo-external forces coming from the disk at the clamped roots of blades at s¼0 in the lateral and longitudinal directions of the respective fan airfoils as, and, The dynamic torque on the shaft at Z ¼ ' due to tangential component of [F Z ] s ¼ 0 is The dynamic torque shown in Eq. (11) causes speed change in the shaft and its rate written as, Additionally in Eqs. (6) and (7), € z 0 and € Z 0 represent the components of clamped-end support accelerations experienced by individual blades in the longitudinal and lateral directions of the blades (see Fig. 2).
The chord length C of the fan-blade airfoil is the distance from its leading to the trailing edge, which at the clamped-end root becomes equal to C ¼ D L =cosb 0 . Thus, in Eqs. (6) and (7) for local support accelerations being fed from the disk to each of the fan blade with its trailing edge located at Z ¼ ', we have, and, The centrifugal stress in the blade at any typical span location 's' from the root of the airfoil is given by, The centrifugal force F cf from each individual blade reacts at the airfoil root, which acts like an external steady-state force in the local span direction on the top of the disk and its magnitude with blade mass M b ¼(r b ChL) is expressed as, The chord-wise component of centrifugal forces F cf on individual blades generate steady-state torque about the span axis of the airfoil due to pre-twist (b L À b 0 ) in the blade, the tendency of which is to untwist the free tip-end of the cantilevered configuration, and is determined as, In the consistent-mass matrix formulation, the pseudo-external forces on the shaft used in Eqs. (2) and (3) from the root of the airfoils at s¼ 0 are, and, However, in the simplified approach of typical rotordynamics applications [17] with lumped-mass matrix formulation, shown above in Eqs. (18) and (19), we use only stiffness terms as, and, The above quasi-external force components on the shaft (F X ) j and (F Y ) j are in fact internal forces in the bladed-disk system, and hence balance out completely by equal and opposite reactions at the clamped-end root reaction of the individual blades. The following Fig. 3 illustrates the asymmetricity in the rotor introduced due to missing blade in one of the disk-sectors. In this figure originally, there are total N blades on the fan rotor in the normal running condition with a fully-balanced rotor-shaft. It is assumed that during a blade-out event one of the blades on the rotor, identified as Blade-N circumferentially located at an angle y 0 , gets suddenly released from the fan-disk.
The torque shown in Eq. (17) caused by the unbalanced centrifugal force F U about the blade local span-axis produces bending moment in the shaft, which balances out completely in an axis-symmetric fully-bladed rotor. However, in a nonsymmetric bladed rotor with a missing bladed, it is not balanced out and acts like an external load in the form of a bending moment at the free overhung end of the shaft. The components of this unbalanced bending moment (BM U ) at any instant due to missing blade during the initial state of (t ¼0) located at y 0 in the global frame of reference are: Due to this imbalance force F U acting radially outward at the heavier side of the disk and simultaneously rotating at the instantaneous angular velocity O of the rotor, the rotor starts moving eccentrically. The local oscillatory movements of the blade root in the longitudinal and lateral directions are given by, and, The objective of the current study is to determine the transient dynamic response of this system, which is subjected to an impulsive load due to sudden blade release and subsequent sustained rotation during windmilling. Now, we will introduce the following Rayleigh-Ritz trial functions, which naturally satisfy all the geometric boundary conditions in the earlier system of Eqs. (2) and (3) for the shaft and (6) and (7) for blades: We define f i ¼ 2iÀ1 ð ÞpÞ=2' À , j i ¼(((2i À1)p)/2L)C 2 , and, a j ¼(((j À 1)p)/C)C 1 with nondimensional terms for the geometrical parameters of the airfoil as C 1 , C 2 (see Appendix B). Then using Eq. (29) and taking its first derivative with respect to s, we also obtain, z, s s,t By substituting the above-mentioned displacement functions and applying the Rayleigh-Ritz technique we obtain the equations of motion in a matrix form made up of typical M, Cand K matrices. In a gyroscopic system, the velocitydependent matrix C is usually broken into two parts called C D for the damping effect and C G for the gyroscopic effect, with former mostly being a symmetric matrix due to damping in the bearings and joints along with the latter as skew-symmetric containing the gyroscopic terms, which results in, Similarly, components of stiffness matrix such as K s contains an elastic modulus-dependent stiffness, whereas K F and K O represent the force-dependent and angular velocity dependent stiffness terms, respectively. The different sub-matrices forming the above global matrix are shown in detail in Appendix B. Here, column vector p(t) represents the external forces which are essentially generated due to friction forces due to blade tip-rub (see Fig. 4).
In the consistent mass matrix formulation, the blade mass terms contributions are weighted according to the deformation shape function Z(s,c,t) in the lateral direction, which results into rigid-body translational inertia mass effect at the airfoil clamped root being as: Similarly in the longitudinal direction of the blade due to deformation function z(s,t), the respective rigid-body translational inertia mass effect at the clamped-end becomes: in addition, the flexural bending of the blades introduces additional inertia terms as, (a) Blade rotary inertia due to tangential motion: (b) Blade gyroscopic inertia due to axial motion: Thus, the global mass and stiffness matrices M and K for the current flexible bladed-rotor system can be expressed in terms of multiple sub-matrices. The global stiffness matrix K is made up of smaller sub-matrices such as, ðZ,U0Þ ,K ðjÞ ðZ,V0Þ and K ðjÞ ðZ,ZÞ etc. Each term of the sub-matrices such as K ðjÞ ðZ,U0Þ ,K ðjÞ ðZ,V0Þ and K ðjÞ ðZ,ZÞ etc. are derived individually and is shown in Appendix B. Using these sub-matrices and the computed values of generalized coordinates like u, v and y (j) , z (j) etc. for any typical blade, the force-displacement relationship in the pseudo-static condition can be written as, In a steady-state condition of a typical rotating blade no. 'j', the external force term p ðjÞ z is due to centrifugal forces and p ðjÞ Z is caused by the blade tip-rub. Using the superscript (j) for the fan blade sequence number (j ¼1yN), the corresponding rotordynamic equations of motion in the matrix form for the flexible bladed-rotor system under consideration can be In the above global matrix, the sign '' is used to indicate the generic form of fully-populated sub-matrices with nonzero terms such as, K ðjÞ BS , K ðjÞ SB , and K ðjÞ B with superscript 'j' ranging from 0 to N for an axi-symmetric rotor, and 0 to (NÀ 1) for nonsymmetric rotor with a blade missing due to blade-loss.
In the above relationships, the stiffness sub-matricessuch as K ðjÞ BS and K ðjÞT SB , due to cross-coupling effect between rotorshaft and the blades, are actually time-dependent with sinusoidal functions containing terms like sin W j þ Ot À Á and cos W j þOt À Á . It should be obvious that the above sub-matrices pertaining to individual fan blades, are repeated along the diagonal of the global stiffness matrix K as many times as the number of fan blades attached to the rotating shaft. An identical process is used to form the respective global inertia matrix M and the velocity-dependent coefficient matrix C. Similarly in a nonsymmetric bladed rotor, we have time-dependent terms appearing in the mass sub-matrix M S with parameters like sin2 W j þ Ot À Á and cos2 W j þ Ot À Á .

Displacement and external force vectors for the rotor-shaft due to blade-loss
In Eq.  (40) The external forces p(t) in such a rotor-shaft are invariably caused due to inherent unbalance in the rotating system. Usually these imbalance forces during the normal operations of the rotor are relatively small in magnitude, and as such do not affect the overall dynamical characteristics of the system. However, sudden loss of a blade from the overhung fan not only generates extremely high impulsive forces on the rotor system, rather it changes the global characteristics of the dynamical system by making the rotor nonsymmetric. In a nonsymmetric bladed-fan rotor, the terms of mass matrix M are no longer constant rather they oscillate with twice the frequency of the spin. Due to large imbalance, the tips of the fan blade also start rubbing against the fan case, which introduce additional set of highly nonlinear contact-impact forces into the system. In this sub-section, we will attempt to account for all these external forces in this typical nonsymmetric rotor system. As external forces on the dynamical system under consideration, suppose we have on a typical jth blade radial interference with the abradable material or the inner surface of the outer casing¼e j (t) such that the radial cold clearance ¼g, which yields, The dynamic rub forces F ðjÞ Rub acting on the blade-tip being contact-impact type in nature will always be compressive with nonzero values only for positive values of e j (t). Additionally, the steady-state centrifugal force F cf on the blade is treated like an external force acting radially outward at the mass-center of the blade with total radial force being, Thus, the respective external force vector terms on the right-hand side of the equations for this dynamical system after losing a fan-blade are:  (43) The tip of a spinning blade located at s ¼L pushes against a rub-surface made up of an abradable material, the radial stiffness of which can be defined as K Rub resulting into the rub forces on the blade tip called F ðjÞ Rub acting along the longitudinal direction of the blade. Thus in the numerical integration scheme, the magnitude of the radial rub-force on the blade tip at the new time-step is determined by the local radial incursion e j (t) using a penalty-method approach. If a typical jth fan blade is rubbing against the abradable on the inside surface of the outer containment ring after consuming the running tip clearance g then the radial component of the rub-force F ðjÞ Rub is, F ðjÞ and, the corresponding tangential component of the rub-force F ðjÞ T will be either, (a) with blade-tip skating/sliding on abradable: (b) with abradable removal across the entire chord at the blade-tip: In the numerical time-marching forward scheme, the blade total tip movement at s¼L in the global frame of reference is written as, In the above equation, the generalized deflection terms (U 0 ,V 0 ) depend upon the bearing support stiffness, whereas the shaft deformation terms (u q ,v q ) depend upon the its elastic stiffness and the generalized terms Z i depend upon the blade stiffness in the longitudinal direction. Thus, it should be obvious that the actual magnitude of rub forces depend upon the component stiffnesses of the entire dynamical system and not just one individual component.
The oscillating axial force P(t) along the longitudinal axis of the shaft due to imbalance and asymmetric mass matrix M creates parametric resonance in the rotating condition. Due to this parametric resonance the rotating shaft may exhibit large amplitude vibrational responses at several sub-harmonic frequencies. The axial force P(t) on the rotor shaft with (NÀ 1) blades remaining on the rotor, is computed as, The axial force on the fan-shaft P(t)shown in Eq. (48) is resisted by the ball bearing on the rotor-shaft assembly. With subscript 'q' varying in the range of (1rqrm) for the number of shaft modes m, the oscillating axial load on the rotating shaft due to unbalance load on the nonaxisymmetric rotor is expressed as, ð Þ, which in the inertial frame of reference is determined as follows: We also consider that the blade-tip rigid body velocity in the tangential direction due to spinning of rotor is À OR, whereas the blade clamped-root tangential velocity due to shaft oscillation is determined by, The dynamic torque in the fan-rotor shaft of the engine T(t) is expressed in the form of a sum of two components, a quasi-steady-state parameter T 0 and a time-dependent parameter of peak magnitude T t as, where, T 0 represents the driving-torque at the turbine-end and is a load-torque at the fan-end of the rotor. It should be noted that T 0 is due to aero-gas load on the fan blades and is proportional to the square of the running speed, which results in, (T 0 ) New ¼(O New /O Old ) 2 (T 0 ) Old . The magnitude of the time-dependent component of the torque T t is controlled by the rub parameters such as the coefficient of friction m, the depth of the radial incursion e j (t) of the blade-tip inside the abradable surface and the radial stiffness K Rub given as, The general solution of Eq. (31) can be assumed in the form of a Fourier series with the generalized coordinate f(t) expressed as, Then the corresponding velocity and acceleration column vectors _ f ðtÞ and € f ðtÞ are written as, and, where, the Fourier coefficient terms such as b 0 , a q , b q etc. are determined from the initial conditions of the dynamical system. The substitution of above terms in the equation of motion and setting the right-hand-side forcing function terms equal to zero, enables us to establish the stability conditions of the rotor in frequency domain [15,16].
In a forced response analysis, once the components of the generalized column vectors _ f and f have been solved, the other rotordynamic parameters such as whirl frequency _ j of the shaft at any axial location can also be determined as, Thus, the shaft whirl frequency in Hz at the axial location of the fan disk is simplified as

Numerical formulation of bearing supports with gap nonlinearity
The current time marching-forward numerical integration scheme is fully capable of handling multiple N c sets of nonlinear bearings used to support the bladed shaft-disk system including bearings with a radial gap. In forming the global stiffness matrix K the contributions of respective bearing stiffness terms S b outlined in Eq. (1) are updated at each time-step 't' and are considered active only when the local (at Z¼z i ) shaft radial deflection d s (z i ) exceeds the pre-specified value of the corresponding radial gap g i of the bearing in that location, such that Using the penalty-method approach, the components of external force vector on the shaft due to nonlinearity and the gap in the bearing supports in the right-hand side of the equation is determined as, where, the shaft radial deformation angle for a particular bearing located at the axial coordinate of Z ¼z i is determined as, 3. Solution of a sample numerical problem As a numerical example, we will consider a hypothetical overhung fan-rotor running at the top speed of 3000 rev/min. The specific geometrical details of the fan-rotor considered here are listed in Table 1.
In the above example data, the effective abradable stiffness K Rub is estimated as, Furthermore, the general test data based upon static coefficient of friction is 0.15, but the dynamic coefficient of friction is usually an order of magnitude less. Hence, in the current transient simulation with the assumption of light-rub resulting in blade-tips simply sliding/skating on the abradable surface, we have used the value of m¼0.01.

Eigenvalue solutions of the flexible-bladed rotor shaft
We have developed a detailed finite element model of the above typical large commercial engine hypothetical turbofan rotor. This bladed-rotor model with the typical values of elastic parameters of various components of the dynamical system such as blade stiffness, abradable stiffness, nonlinear contact stiffness etc. has been analyzed using a commercial nonlinear finite-element code. The frequencies for different components and the blade-disk-rotor assembly are determined in the stationary condition as well as two different speed conditions. In the spinning condition the rotor rotates at 3000 rev/min as the top speed of the shaft and at 1000 rev/min as the windmilling speed of the asymmetric rotor in the damaged condition. The analytically determined first three frequencies and the corresponding mode shapes (1F, 1T, 2T y etc.) for the fan blade are summarized in Fig. 5.
As shown here, the natural frequencies of a stationary blade for the first 3 modes of vibration (1F, 1T and 2T) are 45.89 Hz, 182.78 Hz and 552.50 Hz, respectively. The natural frequencies in the rotating conditions are also determined at two different speeds highlighting the significant effect of stress-stiffening for the first flexural mode (1F), a relatively small effect for the first torsional (1T) mode and a distinct effect of stress-softening for the second torsional (2T) mode of vibration. We have also determined the shaft vibrational mode frequencies of the fully-bladed rotor. When the blades are treated as rigid, the eigenvalue solution yields first two fundamental modes of shaft vibration as 63.8 and 70.11 Hz. The respective mode shape of the rotor is shown in Fig. 6.
However, when the blades are considered flexible, with its representative elastic Young's modulus value, the first shaft bending mode frequency drops by 8.4 percent to 58.45 Hz. As shown in Fig. 7, there is a very clear participation of the blade first flexural mode shapes with the deformation of the shaft.
Interestingly, for this rotor the bending mode of the shaft in the rotating condition at the top speed of 3000 rev/min again moves up to 62.73 Hz, which is very close to the rotor frequency with the rigid blade assumption. Thus, it shows that in case of an axi-symmetric rotor the assumption of blade being rigid is good enough for determining the critical speed of the rotor. However when the same rotor is analyzed with a missing blade making it asymmetric, the shaft mode at 58. 45  It should be noted that the above frequencies still inherently assume that all the bearing supports are fully active. However, in general by the time the windmilling of the asymmetric rotor with a blade missing occurs, the bearing no. 3, closest to the overhung fan-disk may also be damaged or fused due to huge imbalance in the system. The design philosophy behind the failure of bearing ]3 allowed to fail is discussed in quite detail by VanDuyn [47] in his patent disclosure. For this type of asymmetric rotor with a large imbalance, the critical speeds are best established by determining the forced response of the rotor at varying rotational speed of the shaft. Fig. 8 shows the forced harmonic response of such a damaged rotor-shaft by computing its amplitudes of horizontal and vertical motions at the bearing support location separately for each speed (0-3000 rev/min) being changed at the increment of 1 rev/min. The response plot highlights the existence of two resonant speeds for this type of decelerating damaged rotor, with the first one occurring at 2380 rev/min (39.67 Hz) and the second one at 100 rev/min (1.67 Hz). In general a damaged rotor-shaft due to blade-loss invariably decelerates due to aerodynamic stalling of the turbofan. That is the most common observation after a typical blade-loss.

13
A rapidly decelerating rotor-shaft may exhibit its peak vibrational response within a band of speed slightly lower than its resonant speed computed in this manner. However, in a hypothetical scenario of a fan blade-loss event coupled with the simultaneous mid-shaft failure, the turbine-end of the rotor-shaft may accelerate after the blade-out occurs. The original shaft mode frequency (corresponding to 58.45 Hz stationary and 62.73 Hz in rotating condition of 3000 rev/min with 3 undamaged bearing support system), shown in Fig. 7      known as 'decoupled mode', and the lower one at 1.67 Hz can also be easily identified by solving eigenvalue problem of the fan rotor-shaft supported only on two remaining bearings i.e. Bearing-1 and Bearing-2, respectively.
The fan-rotor starts slowing down immediately after the blade-loss as the rotor shaft decelerates and its rotational speed drops from the peak speed of 3000 rev/min to the state of rest. In most of the cases also per FAA requirements for safe operation of the crippled aircraft, this happens in 15 s or less. The resonant speed of 2380 rev/min of the damaged rotor-bearing system exhibits the fan blades momentarily rubbing against the fan-case/abradable during the transient event as the rotor-shaft during the spool-down passes through the 39.67 Hz excitation frequency. As some of the bladetips may experience intermittent rubs, the rotor-shaft response frequency moves up slightly due to additional load-path in the dynamical system. Thus, effectively the fan-rotor response frequency during the transient event of the spool-down may have slight shift due to nonlinearity of the system. These nonlinear contact-rub effects will definitely influence the dynamic response of the damaged rotor during the windmilling phase in a significant way. The dynamic characteristics of such a nonlinear system can only be predicted with any confidence by using a time-marching numerical scheme for integrating the equations of motion. Section 3.2 discusses the transient response of such a rotor.
The lower resonant speed of 100 rev/min may become important in the windmilling of the damaged rotor during the fly-home mission. At this time since the rotor is very lightly damped and the natural frequencies in the two directions are slightly separated, the rotor can wildly oscillate around 100 rev/min and also at multiples of the windmilling speed.

Numerical results of transient response
The transient solution of the dynamical equation is obtained by using the initial conditions at time t ¼0 due to sudden blade-loss and by representing Eq. (31) in the inverse form. In this paper, the above systems of equations of motion are integrated using a sixth-order Runge-Kutta time-marching forward numerical scheme. In this phase of the analysis, it is assumed that a full fan blade gets released at the top speed of 3000 rev/min. The transient analysis is performed for 3 s from the time of blade release under two separate assumptions: (a) Rigid blade formulation. (b) Elastically deformable blade formulation.
The results indicate that under rigid-blade assumption, the rotor responds much quicker than in an identical blade-loss condition of deformed blade formulation. This is because in the early phase of dynamic event, the fan blades offer much bigger resistance due to inertia, but once it starts responding the eventual rotor movement is significantly larger than the rigid-bladed response. The rotor movement may also depend upon the number of blades coming in contact with abradable, which can change every instant.
The transient dynamic response of the rotor is summarized in terms of varying horizontal and vertical movements of the fan-disk by plotting its orbital motion (see Fig. 9). The magnitude of radial excursions uð'Þ and vð'Þ in the horizontal and vertical planes is changing due to rotor being asymmetric. In the orbit plot it can be seen that as the blade-tip pushes against the abradable surface, it bends essentially in the first flexural mode [45], and it loses contact, and then the next blade may come into contact. Although, the flexural bending appears to be the dominant mode of blade deformation during rub, the actual deformed shape at any instant may be a combination of several different modes such as 1F, 1T, 1EB etc. (see Fig. 5). The numerical scheme keeps track of the contact for each individual blade on the fan-rotor by computing radial location of that particular blade-tip using Eq. (58). It is determined at each time-step by the radial coordinate of the blade-tip in the deformed configuration of the shaft. In a flexible bladed-rotor system with each blade having its own lateral and longitudinal degrees-of-freedom (Z,z), the number of blades rubbing can change at each instant. In the numerical integration scheme and the computer code developed for this purpose does a flawless job in keeping track of all the blades which might potentially rub. In our current example, sometimes only one blade has been found to be rubbing and sometimes at the most two to three blades also. It is worth noting that we do not make any assumption a priori about the number of blades which may or may not rub.
As shown in Fig. 9, the transient nondimensional radial movement of the fan-disk is as much as 5 percent of fan radius from the engine axis. However, under a realistic condition of flexible blade formulation most of the eccentric movement of the fan-disk center is within 72 percent range. The corresponding blade bending movement is about 2 percent, indicating the magnitude of the radial load on the bearing would have peaked at 11.3 M N. Here, the stiffness and damping of the bearing no. 3 is removed out of the transient analysis after the load on the bearing support has reached to its maximum load carrying capacity of 1 M N. Fig. 10 highlights the differences in the transient behavior of deformable blade model versus rigid-blade response of the rotor. Under the rigid blade solution, there is a phase change of the radial excursion of the shaft with the heavier-side of the disk rubbing against the case in the beginning and then changing to the light-side of the disk making contact with the fan-casing. Due to the whirling motion of the shaft the phase-change between the heavy-side of the disk to the light-side of the disk rubbing against the case is an intermittent event. With this back-and-forth phase-change there may be certain instants occurring periodically when for small durations no blades will be rubbing at all. This duration of no contact can be longer if the particular disk-sector with the missing blade is moving closer to the fan casing. In the current example, such a dynamic behavior of the rotor has been exhibited through a purely analytical solution using the sixth-order Runge-Kutta numerical scheme, but this has also been observed in high-speed movie of a typical fan-blade-out event. Under the rigidblade assumption, the rotor moves radially outward much faster, and the computed radial deflection of the shaft peaks at inner surface of the hard containment casing. The dynamic response of the rotor with rigid-blades is very much identical to the results shown by Lahriri et al. [42]. However, the deformable blade solution indicates that blade tip-bending against abradable surface plays a very important role in determining the orbit as well as the rate at which rotor moves off-center. The current pre-twisted blade model automatically takes into account of many different types of complex mode shapes as well (see Fig. 5).

Windmilling of an asymmetric rotor after blade-loss
A typical turbofan rotor starts behaving like a turbine rotor, when it is subjected to a head-wind, which feeds energy to the dynamical system. The safety objective of FAA regulation 33.74 is to ensure that an engine that continues to rotate after shutdown will not create a hazard to the aircraft. Continued rotation refers to a condition in which any part of the main rotating system in an engine continues to rotate after the engine has been shut down. Continued rotation can be caused by windmilling or mechanical effects, or a combination of both. Windmilling is the rotation of a nonoperating engine due to the airflow-induced forces on the blades caused by the forward motion of the aircraft. Mechanical effects include, for example, drive shaft clutch drag in some multiengine rotorcraft installations, which may result in continued rotation of the engine after it has been shut down. The FAA document [8] describes, '' The vibratory loads resulting from the failure of a fan blade have traditionally been regarded as insignificant relative to other portions of the design load spectrum for the airplane. However, the progression to larger fan diameters and fewer blades with larger chords has changed the significance of engine structural failures that result in an imbalanced rotating assembly. This condition is further exacerbated by the fact that fans will continue to windmill in the imbalance condition following engine shut down. Current rules require provisions to stop the windmilling rotor where continued rotation could jeopardize the safety of the airplane. However, large high bypass ratio fans are practically impossible to stop in flight.'' In a damaged fan rotor after losing a blade, intermittent rub can produce very large vibrational response, which may contain many different types of sub-harmonic characteristics. In an aircraft engine after a fan blade-out event, the windmilling of an asymmetric rotor takes place during the fly-home period. A typical turbofan rotor usually will windmill at 1/3 of the full-thrust speed of the rotor. Using the same transient response technique described in Section 3.2, we have analyzed the dynamic response of the damaged rotor, if it were windmilling at 1000 rev/min. The time-domain dynamic response of the damaged nonaxisymmetric rotor at 1000 rev/min is summarized as follows: The above time-domain results are analyzed at 1000 rev/min for 8 s. Here, the forcing function is the unbalanced force caused by the missing blade as well as the solution scheme updates the mass and stiffness matrices as the rotor spins. As shown in Fig. 11, the peak horizontal movement of the rotor at the center of the fan disk oscillates between 0.5 cm and 2.0 cm. The zoomed-in dynamic response of this rotor, shown in Fig. 12, exhibits many different types of frequency content. This time-domain response includes the effect of both oscillations in the mass matrix at twice the rotational speed of the rotor as well as the time-dependent terms in the stiffness matrix due to intermittent blade tip-rub against abradable inner surface.
Similar time-domain response of the rotor in the vertical direction has been plotted and is shown in Fig. 13. The asymmetricity of the rotor is quite evident by the fact that the horizontal direction response is significantly different than its vertical direction response. As shown in Fig. 13, the peak vertical movement of the rotor at the center of the fan disk oscillates between 1.0 cm and 1.7 cm. Fig. 14 shows the zoomed-in view of the vertical response in a quasi-steady state of the rotor windmilling at 1000 rev/min. The time-domain response of the asymmetric rotor during windmilling shown in Fig. 14 clearly highlights the fact that in its dynamic behavior, there are multiple frequencies superimposed on one another. We also recover the rub response of the running asymmetric rotor during windmilling in the aftermath of the blade-loss. At this point, we will determine the frequency content of the dynamic response by using a fast-Fourier transform of the time-domain data displacement data. The actual time-domain behavior of the rotor-shaft comprises of responses at several different frequencies, and the instantaneous whirl-frequency of this type of asymmetric fan-rotor can also be determined using Eq. (58).
The above time-domain results also help us in determining the rub-load on the rotor coming from the blade-tips rubbing against the abradable. Fig. 15 shows the time-domain instantaneous magnitudes of the analytically determined rub-loads on the fan-rotor. The plot shows that the magnitude of the contact-impact rub-load in the radial direction of this rotor at the fan-disk oscillates about the mean value of 100,000 N with its amplitude peaking up to 190,000 N. A simple hand-calculation shows that an impulsive load of this magnitude in the longitudinal direction of the airfoil will produce a compressive membrane strain in the vicinity of rub-zone at the blade tip of the order of 0.033 percent. This strain value is very well within the range of measured strain data reported earlier [11,46]. The detailed comparison of the current analytical results versus measured data is discussed later on in Section 5 of this paper. Fig. 16 shows that during windmilling a damaged nonsymmetric rotor would respond at a sub-harmonic frequency. For example, the natural frequencies of this rotor at 1000 rev/min is at 70.35 Hz and 62.85 Hz, but the rotor shows high magnitude oscillation not only at its excitation frequency of 16.66 Hz, rather it also exhibits sub-harmonic response at 8.25 Hz. The sub-harmonic frequency coincides with half of the excitation frequency O/2, and the system responds at its sum-frequency of 25.4 Hz also.    Fig. 17 is plotted for a steady-state rotational speed of 1000 rev/min (16.67 Hz) as a result of unbalance caused by the missing blade. Any whirl response of the rotor at a frequency higher than the running excitation frequency will result in a forward-whirl response of the rotor. Similarly, the dynamic response occurring at less than the running excitation frequency will exhibit backward whirl characteristics of the rotor shaft. The forward and backward whirl response of the fan-rotor can very easily be visualized by creating an animation of the fan-disk orbit from the transient response plotted in Figs. 11 and 13. In addition, Fig. 17 shows that this damaged asymmetric rotor during windmilling at 1000 rev/min will go into distinct backward-whirl behavior at 6.1 Hz and other two minor ones at 10.9 Hz and 12.2 Hz. In addition, it will show a complex forward-whirl behavior with its whirling frequencies at 27.9 Hz and 34.0 Hz. There are several other minor blips in the whirl frequency plot indicating that this shaft may exhibit a very complex forward-whirl behavior at other frequencies such as 21.6 Hz, 24.1 Hz, 38.3 Hz and 44.6 Hz etc. as well. Depending upon the magnitude of the shaft oscillation during whirling, the radial rub-load on the fan rotor will be oscillating between 25,000 and 190,000 N (see Fig. 15). The most dominant whirl characteristic for this long rotor will be a synchronous whirl in the vicinity of excitation frequency of 17 Hz. The existence of multiple whirl frequencies indicates that the different sections of the fan-rotor will be whirling at varying frequencies. Considering its critical speed at 100 rev/min (see Fig. 8), the low frequency whirl may result into large-magnitude oscillations of the shaft during windmilling.   This rotor depending upon the blade tip clearance g will have occasional rub in some isolated circumferential locations of the fan disk. During the solution of the governing equations numerically, we keep track of the particular blade tips' coordinate locations to determine whether any parts of the blade tips are in contact with the abradable on the containment casing. The time-depending circumferential location of a typical blade identified as 'Blade-j' which might be rubbing at any time 't' is expressed as, (W j þOt). The intermittent rubbing with blades being flexible in the longitudinal direction introduces a coupling between the longitudinal and transverse deformations along its span. This coupling is caused by the Coriolis forces as a result of longitudinal wave of the blade traveling up-and-down along the span due to the blade rotating in the centrifugal force field [46].

Comparison of measured rub data with analytically computed results
In order to establish the accuracy of the current analytical model for its time-domain results, the numerically computed transient dynamic response of a typical high-pressure compressor blade is compared with the previously published [11,46] rig-test measured data under controlled conditions. In the current numerical scheme, the transient analysis has been carried out for 4 repeated impacts with 0.1 mm radial incursion with the blade tip tangential velocity being at 400 m/ s. The measured dynamic data from a spanwise strain gage at the fillet of the airfoil root is compared with the numerically computed strain time-history near the clamped-end of the airfoil. From the two sets of plotted data illustrated in Fig. 18, it can be seen that the transient analytical results predicts the dynamic characteristics and the resulting strain time-history in the rubbing blade very well.

Conclusions
The present study involves a comprehensive theoretical methodology to simulate a very complex nonlinear dynamic event. The goal is to address the theoretical complexities and develop the proper mathematical basis so that such an analysis can be carried out accurately and the results can be evaluated with full confidence. The analytical method is properly complemented by a detailed numerical simulation.
In this paper, we have derived all the rotordynamic equations of a nonsymmetric rotor with deformable blades, which can rub intermittently as a post-blade-out with the outer casing. The paper presents an analytical method to illustrate the nonlinear dynamic effect of blades rubbing against the rigid outer case in a rotating machinery. The equations are  represented in a classical matrix form of M, C, and K matrices. Each and every term of these matrices have been derived independently and have been shown in the appendix of the paper. Both mode shapes as well as transient response of the rotor is determined over a complete spectrum of the speed-range. An attempt has been made to quantify the magnitude of the contact load during windmilling of a damaged rotor with a blade missing, and it is shown that during the hard rub against the outer case, the sudden contact-impact load can go up by an order of magnitude over the unbalance load. For a typical rotor blade the radial load from the case to the blade tip can last for 0.05 m s with a peak magnitude reaching up to 0.1 M N. However, it is recognized that the analytically predicted dynamic contact loads in this paper provide an upper bound of the axial load magnitude on the blades due to following two underlying assumptions: (1) Plasticity is not considered and blades remain elastic for the entire duration of impact.
(2) Geometrically blades are considered like a cantilever beam of uniform cross-section in the undeformed configuration.
Blades are allowed to deform in the lateral direction either due to tangential friction load at the tip, or pulse-buckled deformation under axial impact load. (3) In the dynamic response of a pre-twisted blade subjected to contact-impact rub-loads, the fan-blades can deform in several different complex modes, such as edge-wise bending, if the contact takes place only at extreme chordlocations, such as at the leading or trailing edge of the blade.
Due to these underlying assumptions, the actual contact load is expected to be somewhat lower than computed here for elastic condition. The method can be refined to take into account of the yield stress and elastic-plastic deformation of the blade material.
where the time-dependent gyroscopic and inertia terms contributed by each individual blade (j ¼1yN b ) on the fan-disk rotating about the engine axis result into its cumulative sum as, (a) Bladed-disk dynamic moment at the shaft center about the global X-axis: and, (b) Bladed-disk dynamic moment at the shaft center about the global Y-axis: Appendix B. Determination of terms in conventional M, C and K matrices For brevity, we introduce coefficients such as a i , j i , f i etc. to represent the sinusoidal terms of the series to develop the deformation shape functions for various components of the turbofan rotor used in the current analysis. Here, we define, where nondimensional terms are, Using the nomenclatures shown above, we have the Rayleigh-Ritz trial functions described in Eqs. (26)- (29) for the shaft and the blade assumed as, (a) Shaft deformation (u,v): with (0rs rL) and, ( À C/2rc rC/2), where In Eq. (B3), the coefficients c 2 and c 3 are dependent upon the blade boundary conditions. Thus, on applying the boundary conditions of bending moment and shear forces at the free-end of the cantilever blade being zero with W 00 i ðLÞ ¼ 0 and, W 000 . Thus, the term W i (s) for the lateral deflection of the blade with all boundary conditions satisfied is, Hence, the transient transverse deformation term Z¼Z(s,c,t) of any typical blade for any point on its mid-surface located at (s,c) in its local shell-coordinates, is obtained as, Here, for capturing all relevant blade vibrational modes we have considered n1 number of terms in the span direction and 4 terms of the series in the chord or width direction [44] of the blade. Similarly, the transient longitudinal deformation component z¼z(s,t) of any typical blade is expressed with n 2 number of terms of the series as, The typical nonzero terms of the global M, C and K matrices are shown in Eqs. (B8)-(B64). In these equations, the superscript '(j)' represents the fan blade sequence number with j¼ 1 being the first trail blade on the fan rotor after the blade-loss. In an attempt to keep the dynamical sub-matrices derivations general enough for any rotating blade positioning we have expressed respective vectors as sin(W j þOt) and cos(W j þOt) components. However for specific turbofan blade-loss applications, we can write Ot¼0, which enables us to simplify the corresponding vectors in terms of sin W j and cos W j components, respectively. The highest limit on the value of j can be 'N' as the number of blades on the fully-bladed fan rotor-shaft, and (N À1) being as the as the number of blades left on the damaged fan rotor. In the following equations, N b is the number of active blades on the fan-disk, which is usually an even number N in a regular symmetric rotor. In the aftermath of a blade-loss with (NÀ 1) remaining blades on the rotor-shaft, the fan-rotor becomes nonsymmetric. It should be noted that in the dynamical equations for the shaft-disk-bladed-rotor system the number of rows and columns for various sub-matrices depend upon the following numbers: m ¼number of deflection shape function terms used to define shaft deformation (u,v) n 1 ¼number of deflection shape function terms used to define blade lateral (Z-degree of freedom) deformation and n 2 ¼number of deflection shape function terms used to define blade longitudinal (z-degree of freedom) deformation.
In addition, subscripts (p,q) refer to the terms in pth row and qth column for respective sub-matrices, shown below.
(B1) Inertia terms of acceleration-dependent sub-matrices M p,q The contribution of inertia effect of the fan-blade movement towards the rigid-body shaft motion coming through the airfoil clamped-end can be accounted as an equal weighted -average of the number of terms used in the blade deformation (see Eqs. (28) and (29)). Thus, the following sub-matrices shown in Eqs. (B10) and (B11) used in the consistent mass matrix formulation will be repeated j ¼1yN b times for each blade: In the following Eq. (B17), N b is the number of active blades on the fan-disk, which is usually an even number N in a regular symmetric rotor. In the aftermath of a blade-loss with (NÀ 1) remaining blades the fan-rotor becomes asymmetric, which makes the mass-matrix term [M (U,U) ] p,q time-dependent resulting into twice of the running speed frequency excitation due to contribution of the cos 2 W j þ Ot À Á term.
The following 2 sub-matrices of Eqs. (B19) and (B20) are set to zero under lumped-mass assumptions but can be used in the consistent-mass matrix formulation for each individual airfoil by varying j ¼number of fan blades (N b ) times:  (B21) In the following Eq. (B22), N b is the number of active blades on the fan-disk, which changes its value for the reasons described earlier following Eq. (B16): The following set of sub-matrices of Eqs. (B23) and (B24) are set to zero under lumped-mass assumptions but can be used in the consistent-mass matrix formulation for each individual airfoil by varying j¼number of fan blades (N b ) times: (B2) Damping and circulatory terms of velocity-dependent sub-matrices C p,q The velocity-dependent coefficient matrix C includes the effect of gyroscopic and viscous damping present in the dynamical system and is expressed as, The components of gyroscopic matrix term C G are always skew-symmetric in nature and contain coefficients like 72OM to account for the spinning rotor with,