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

Sunil Sinha

► To cite this version:

Sunil Sinha. Rotordynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads. Journal of Sound and Vibration, 2013, 332 (9), pp.2253-2283. 10.1016/j.jsv.2012.11.033. hal-01555281

HAL Id: hal-01555281

https://hal.science/hal-01555281

Submitted on 7 Jul 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Rotodynamic analysis of asymmetric turbofan rotor due to fan blade-loss event with contact-impact rub loads

Sunil K.Sinha

GE Aviation, General Electric Company, Cincinnati, OH 45215, USA

Loss of a blade from a running turbofan rotor introduces not only huge imbalance into the dynamical system rather it makes the entire rotor asymmetric as well. In a nonsymmetric rotor, the various terms of mass, gyroscopic and stiffness matrices also become time-dependent. In this paper, all the dynamical equations include the effect of the rotary inertia and gyroscopic moments as a result of both shaft bending as well as staggered blades flexing in-and-out of the plane of the disk. The governing equations also account for internal material damping in the shaft and the external damping in the support bearing system. In addition to the unbalance load at the disk location, the shaft may also be subjected to a torque and axial forces. Here, the fan blades are modeled as pre-twisted thin shallow shells. They have coupled flexural-torsional motion in the lateral out-of-plane direction as well as extensional degrees-of-freedom in the longitudinal spanwise direction of the blade airfoil. The effect of blade tip rub forces being transmitted to the shaft are analyzed in terms of the dynamic stability of the rotor, especially during windmilling.

## 1. 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-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
blade-loss event on an engine are of proprietary nature and are not currently available in public domain. Any experiment simulating blade-loss or tip-rub under controlled conditions [11] is not only extremely expensive, rather sometimes they involve serious risk to the safety of the rig and the test personnel as well. As such there is a scarcity of test data, to compare the results of any fan blade-out analytical model in the published literature. However there are ample information from the actual field events and their investigation reports [12], which provide the necessary evidence about the magnitude of the dynamic loads due to fan blade-loss under operating conditions of a turbofan engine.

In the aftermath of blade-loss, the fan-disk rotating system starts moving off-center in the heavier side of the disk due to sudden large imbalance introduced on the rotor-shaft. This eccentric rotor no longer remains symmetric about its center of rotation, and the spinning of the asymmetric rotor introduces another set of nonlinearities. Recently an attempt has been made to analyze the free vibrational characteristics of a weakly asymmetric rotor [13]. The present work is the first paper of its kind in the field of rotormechanics which addresses the forced transient response of an asymmetric flexible bladed-rotor with also considering the casing rub-surface flexibility in its formulation while the outer fan-case is still effectively undeformed. The inner surface abradable flexibility is expressed by an equivalent term $K_{\text{Rub}}$ and its dynamic effect on the blade-tip is included using a penalty method approach.

In this paper, we consider the transient dynamic response of a complete rotor system in a jet engine with its blade rubbing against the containment case. In the past, researchers have studied the dynamic stability of a typical rotor shaft in many different forms by varying parameters such as internal damping, torque, tip-rub, pulsating axial force etc. [14–25]. However in general, these studies invariably consider the bladed-disk assembly as rigid with entire flexibility of the dynamical system attributed to the shaft material with an artificial nonlinear spring to account for blade-tip rub forces. Flexibility of a blade-disk assembly has always been a concern in terms of dynamic behavior of a bladed-disk-rotor system, but no comprehensive studies have been carried out which has focused on all the interaction effects. The blade/casing interaction can be of many different types, such as deformable case response [26] or a plastically deformable abradable inner surface [27,28] with material removal being considered due to tip-rub. In 1995 Hsieh and Abel [29] were the first ones to use finite-element analysis to investigate into the combined effect of flexibilities of disk and individual blades. Genta [30] carried out a stability analysis of a bladed-rotor dynamical system. For the first time, Sinha [31] combined the bending flexibility of individual blades with that of the rotating shaft in a fixed coordinate system by treating each component as separate Euler-Bernoulli's beams. Lesaffre et al. [32] extended Sinha's model to carry out the dynamic stability analysis of a typical overhung disk-rotor system commonly used in a turbofan engine. Similar approach was followed by Cao et al. [33] to analyze the nonlinear effect of such a system with a squeeze film damper. Legrand and his team's [34,35] research have been focused on the interaction effects of flexibility of the bladed-disk and the outer casing. Gruin et al. [36] studied the nonlinear dynamics of a bladed dual shaft. Recently, Shang et al. [37] have used a simplified rotor-stator model for analyzing the global response during the rub with friction. Arzina [38] wrote a thesis analyzing the compressor blades rubbing against the case. Williams [39] used an implicit time marching forward numerical scheme to analyze heavy tip-rub. Turner et al. [40] considered the blade deflection characteristics during the rub event. Torkhani et al. [41] experimentally and numerically analyzed different amounts of rub during speed transients of the rotor. Lahriri et al. [42] studied the stator-rotor contact dynamics using a nonideal drive system. Kim and Lee [43] have done a groundbreaking work by studying the interaction effects between an asymmetric rotor with an anisotropic stator.

Here in the current formulation, the blade geometry is assumed to be a pre-twisted airfoil [44] which can deform both in the flexural as well as in extensional modes [45]. Both the shaft and all the individual blade deformations are represented by a set of Rayleigh-Ritz's trial functions which satisfy all the geometric as well as natural boundary conditions arbitrarily. Thus, by using Galerkin's method a direct matrix formulation of the nonlinear dynamic problem has been developed, which does not make any such simplifying assumptions, nor does it need any transformation. Yet, it is general enough to include the effect of viscoelasticity in the shaft material, aerodynamic load on the blades, nonuniform support damping as well as the gyroscopic terms. In the prevailing rotordynamics approach, due to complexity of the nonlinear analysis, the flexible shaft-disk of the dynamical system are frequently reduced to much simplified two-degrees of freedom Jeffcott-rotor models made up of a set of lumped-masses and springs. Furthermore, due to large computation time involved, attempt is usually limited to obtain the solution only in frequency domain, where peak response are determined at certain pre-determined harmonic excitation frequencies, which happens to be the entire range of the running speed of the rotor during normal operation.

## 2. 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
<figure><img src="image_2.png" /><figcaption>Fig. 1. Overhung fan disk rotor shaft being supported by multiple bearings with a typical pre-twisted airfoil fan blade.</figcaption></figure>

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 $\mu F_{Rub}$ acting in the direction opposite to the direction of shaft rotation. The effective coefficient of friction $\mu$ 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 $l$, which may be either ball-bearing or roller-bearing, with their respective stiffness $S_b^{(i)}$ and damping $D_b^{(i)}$ 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,

$$ S_b^{(i)} = \begin{bmatrix} S_{xx}^{b(i)} & S_{xy}^{b(i)} \\ S_{yx}^{b(i)} & S_{yy}^{b(i)} \end{bmatrix}_{\substack{B_{xx}^{b(i)} \\ B_{yy}^{b(i)}}} \quad \text{and,} \quad D_b^{(i)} = \begin{bmatrix} D_{xx}^{b(i)} & D_{xy}^{b(i)} \\ D_{yx}^{b(i)} & D_{yy}^{b(i)} \end{bmatrix} \quad (1) $$

where, matrix terms such as $S_{xx}^{b(i)}$, $S_{yy}^{b(i)}$, $S_{xy}^{b(i)}$, and $S_{yx}^{b(i)}$ have the dimensions of linear stiffness as force per unit length, and $B_{xx}^{b(i)}$, and $B_{yy}^{b(i)}$ have the dimensions of bending stiffness as moment/unit radian. For, roller bearings with no rotational resistance the terms $B_{xx}^{b(i)}$ and $B_{yy}^{b(i)}$ will be zero, but for ball bearings with rotational load carrying capability the terms $B_{xx}^{b(i)}$ and $B_{yy}^{b(i)}$ 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 = \ell$. 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:

$$ 
\begin{aligned}
& \left[ E_s I + 2\psi H \frac{\partial}{\partial t} \right] u,_{ZZZZ} + \rho_s A_s \ddot{u} - I_d \ddot{u},_{ZZ} - 2I_d \Omega [\dot{v},_{ZZ} - \Omega u,_{ZZ}] + 2\psi H \Omega v,_{ZZZZ} + T(t)v,_{ZZZ} \\
& -P(t)u,_{ZZ} + \sum_{i=1}^{N_c} [S_{xx}^{b(i)} u + S_{xy}^{b(i)} v + B_{xx}^{b(i)} u,_{ZZ} + D_{xx}^{b(i)} \dot{u} + D_{xy}^{b(i)} \dot{v}] \delta(Z-Z_i) \\
& + [(M_D)_0 \ddot{u} - (J_{D-YY})_0 \ddot{u},_{ZZ} - (J_{D-XX} + J_{D-YY})_0 \Omega [\dot{v},_{ZZ} - \Omega u,_{ZZ}]] \delta(Z-\ell) \\
& -\Omega^2 \left[ (\rho_s A_s) u + \sum_{k=1}^{N_D} (M_D)_k u \delta(Z-Z_k) + (M_D)_0 u \delta(Z-\ell) \right] \\
& -\sum_{k=1}^{N_D} 2 (J_D)_k \Omega [\dot{v},_{ZZ} - \Omega u,_{ZZ}] \delta(Z-Z_k) \\
& + \sum_{j=1}^{N_b} M_b \left[ \ddot{u} - (r+L/2)^2 [\ddot{u},_{ZZ} + 2\Omega \dot{v},_{ZZ} - \Omega^2 u,_{ZZ}] \cos^2(\vartheta_j + \Omega t) \right. \\
& \left. - (r+L/2)^2 [\ddot{v},_{ZZ} - 2\Omega \dot{u},_{ZZ} - \Omega^2 v,_{ZZ}] \sin(\vartheta_j + \Omega t) \cos(\vartheta_j + \Omega t) \right] \delta(Z-\ell) \\
& + \sum_{j=1}^{N_b} \left[ -CD\eta,_{sss} \cos \beta_0 \sin(\vartheta_j + \Omega t) + [\eta,_{ss} \cos \beta_0 - r\eta,_{sss} \sin \beta_0] CDu,_{ZZ} \cos(\vartheta_j + \Omega t) \right] \\
& \qquad \left. \delta(Z-\ell) \right]_{s=0} \\
& = [F_U \cos(\Omega t + \pi + \theta_0)] \delta(Z-\ell - D_L/2)
\end{aligned}
\quad (2)
$$

$$ 
\begin{aligned}
& \left[ E_s I + 2\psi H \frac{\partial}{\partial t} \right] v,_{ZZZZ} + \rho_s A_s \ddot{v} - I_d \ddot{v},_{ZZ} + 2I_d \Omega [\dot{u},_{ZZ} + \Omega v,_{ZZ}] - 2\psi H \Omega u,_{ZZZZ} - T(t)u,_{ZZZZ} \\
& -P(t)v,_{ZZ} + \sum_{i=1}^{N_c} [S_{yx}^{b(i)} u + S_{yy}^{b(i)} v + B_{yx}^{b(i)} v,_{ZZ} + D_{yx}^{b(i)} \dot{u} + D_{yx}^{b(i)} \dot{v}] \delta(Z-Z_i) \\
& + [(M_D)_0 \ddot{v} - (J_{D-XX})_0 \ddot{v},_{ZZ} + (J_{D-XX} + J_{D-YY})_0 \Omega [\dot{u},_{ZZ} + \Omega v,_{ZZ}]] \delta(Z-\ell) \\
& -\Omega^2 \left[ (\rho_s A_s) v + \sum_{k=1}^{N_D} (M_D)_k v \delta(Z-Z_k) + (M_D)_0 v \delta(Z-\ell) \right] \\
& + \sum_{k=1}^{N_D} 2 (J_D)_k \Omega [\dot{u},_{ZZ} + \Omega v,_{ZZ}] \delta(Z-Z_k) \\
& + \sum_{j=1}^{N_b} M_b \left[ \ddot{v} - (r+L/2)^2 [\ddot{v},_{ZZ} - 2\Omega \dot{u},_{ZZ} - \Omega^2 v,_{ZZ}] \sin^2(\vartheta_j + \Omega t) \right. \\
& \left. - (r+L/2)^2 [\ddot{u},_{ZZ} + 2\Omega \dot{v},_{ZZ} - \Omega^2 u,_{ZZ}] \sin(\vartheta_j + \Omega t) \cos(\vartheta_j + \Omega t) \right] \delta(Z-\ell) \\
& + \sum_{j=1}^{N_b} \left[ -CD\eta,_{sss} \cos \beta_0 \cos(\vartheta_j + \Omega t) + [\eta,_{ss} \cos \beta_0 - r\eta,_{sss} \sin \beta_0] CDv,_{ZZ} \sin(\vartheta_j + \Omega t) \right] \\
& \qquad \left. \delta(Z-\ell) \right]_{s=0} \\
& = [F_U \sin(\Omega t + \pi + \theta_0)] \delta(Z-\ell - D_L/2)
\end{aligned}
\quad (3)
$$

where, the time-dependent gyroscopic and inertia terms on the shaft are the cumulative sum of the blade sequence number ranging from $j=1...N_b$ attached to the fan-disk and rotating about the engine axis. Thus $N_b$ makes a step-change with $N_b=N$ for a fully-bladed normal running rotor with its value suddenly changing to $N_b=(N-1)$ for a damaged turbofan-rotor after losing a blade from the disk. Similarly:

(a) Bladed-disk dynamic moment at the shaft center about the global X-axis:

$$ M_X = \sum_{j=1}^{N_b} M_b (r+L/2)^2 \left[ \begin{aligned} & (\ddot{v},_Z - 2\Omega \dot{u},_Z - \Omega^2 v,)_Z \sin^2(\vartheta_j + \Omega t) \\ & + (\ddot{u},_Z + 2\Omega \dot{v},_Z - \Omega^2 u,)_Z \sin(\vartheta_j + \Omega t) \cos(\vartheta_j + \Omega t) \end{aligned} \right]_{Z=\ell} \quad (4) $$

and,

(b) Bladed-disk dynamic moment at the shaft center about the global Y-axis:

$$ M_Y = \sum_{j=1}^{N_b} M_b (r+L/2)^2 \left[ \begin{aligned} & (\ddot{u},_Z + 2\Omega \dot{v},_Z - \Omega^2 u,)_Z \cos^2(\vartheta_j + \Omega t) \\ & + (\ddot{v},_Z - 2\Omega \dot{u},_Z - \Omega^2 v,)_Z \sin(\vartheta_j + \Omega t) \sin(\vartheta_j + \Omega t) \end{aligned} \right]_{Z=\ell} \quad (5) $$

In the above Eqs. (2) and (3), $T(t)$ is the dynamic torque in the shaft, which has positive sign if it is acting in the direction of the spin of the shaft and would have negative sign otherwise. Thus, for equations written earlier for the shaft with
external forces coming from the clamped roots of individual blades in the lateral and longitudinal directions for each blades are, $-[F_\eta]_{s=0}$ and $-[F_\zeta]_{s=0}$, respectively. Here again, the blade-count parameter $j$ varies from 1 to $N$ in a fully-bladed 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 $(\beta_L - \beta_0)$ radians with a negative sign. It should be noted that the sample blade shown in Fig. 1 has a positive $\beta$, 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 $\eta(s,c,t)$ in the thickness direction and can also have membrane stretching due to longitudinal deformation $\zeta(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:

$$ 
\begin{aligned}
& D\nabla^4\eta + D(\beta')^2[3\eta_{,ss} + \eta_{,cc}] - D(\beta')^4 \left[s^2\eta_{,ss} + \frac{(1-v_b)}{2}s^2\eta_{,cc} + 2s\eta_{,s}\right] \\
& \qquad - \frac{F_{\text{Rub}}}{C}[\eta_{,ss} - (\beta')^2\eta] - h\frac{\partial(\sigma_{cf})_s}{\partial s}\eta_{,s} - h(\sigma_{cf})_s\eta_{,ss} \\
& \qquad - 2\rho_b h \Omega \cos(\beta_0 + \beta's)\dot{\zeta} - \frac{\mu F_{\text{Rub}} \cos \beta_L}{C}\delta(s-L) - \rho_b h \Omega^2 \cos^2(\beta_0 + \beta's)\eta \\
& \qquad + \rho_b h \ddot{\eta} - \rho_b (h^3/12)(\ddot{\eta}_{,ss} + \ddot{\eta}_{,cc}) + \rho_b h (\beta's)^2[\ddot{\eta} - v_b(Cc/L)^2\ddot{\eta}_{,cc}]_0 = (1/C)[-F_\eta + M_b\ddot{\eta}_0]_{s=0}\delta(s-0)
\end{aligned}
\quad (6) $$

and,

$$ -\frac{E_b}{(1-v_b^2)} \frac{\partial^2 \zeta}{\partial s^2} + \rho_b \ddot{\zeta} + 2\rho_b \Omega \cos(\beta_0 + \beta's)\dot{\eta} - \rho_b [r+s+\zeta]\Omega^2 - [F_{\text{Rub}}/Ch]\delta(s-L) = (1/Ch)[-F_{\zeta} + M_b \ddot{\zeta}_0]_{s=0}\delta(s-0) \quad (7) $$

In the above Eqs. (6) and (7), the force terms $[F_\eta]_{s=0}$ and $[F_\zeta]_{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,

$$ [BM_\eta]_{s=0} = -CD\eta_{,ss} = -\frac{CE_b h^3}{12(1-v_b^2)}\eta_{,ss} \quad (8) $$

$$ [F_\eta]_s = CD\eta_{,sss} = \frac{CE_b h^3}{12(1-v_b^2)}\eta_{,sss} \quad (9) $$

and,

$$ [F_\zeta]_s = 0 = -\frac{CE_b h}{(1-v_b^2)}\zeta_{,s} \quad (10) $$

The dynamic torque on the shaft at $Z = \ell$ due to tangential component of $[F_\eta]_{s=0}$ is

$$ T(t) = r(CD\eta_{,sss})\cos\beta_0 \quad (11) $$

The dynamic torque shown in Eq. (11) causes speed change in the shaft and its rate written as,

$$ \frac{d\Omega}{dt} = \frac{r(CD\eta_{,sss})\cos\beta_0}{(J_{D-XX} + J_{D-YY})} \quad (12) $$

Additionally in Eqs. (6) and (7), $\ddot{\zeta}_0$ and $\ddot{\eta}_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/\cos\beta_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 = \ell$, we have,

$$ \ddot{\zeta}_0 = \left[ \left( \dot{u} + \ddot{u}, z \frac{D_L}{2} \right) + 2 \left( \dot{v} + \ddot{v}, z \frac{D_L}{2} \right) \Omega - \left( u + u, z \frac{D_L}{2} \right) \Omega^2 \right]_{Z=\ell} \cos(\Omega t + \vartheta_j) \\ + \left[ \left( \ddot{v} + \ddot{v}, z \frac{D_L}{2} \right) - 2 \left( \dot{u} + \ddot{u}, z \frac{D_L}{2} \right) \Omega - \left( v + v, z \frac{D_L}{2} \right) \Omega^2 \right]_{Z=\ell} \sin(\Omega t + \vartheta_j) \quad (13) $$

and,

$$ \ddot{\eta}_0 = \left[ \begin{aligned} & [\ddot{u}, z + 2\dot{v}, z \Omega - u, z \Omega^2] \sin\beta_0 \\ & - \left[ \left( \ddot{v} + \ddot{v}, z \frac{D_L}{2} \right) - 2 \left( \dot{u} + \ddot{u}, z \frac{D_L}{2} \right) \Omega - \left( v + v, z \frac{D_L}{2} \right) \Omega^2 \right] \cos\beta_0 \end{aligned} \right]_{Z=\ell} \cos(\Omega t + \vartheta_j) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic diagram for the fan blade vibratory motion ($\eta_0$, $\zeta_0$) at the airfoil root caused by the shaft and fan-disk movement.</figcaption></figure>

$$ + \left[ \begin{aligned} & [\ddot{v},z-2\dot{u},z\Omega-\dot{v},z\Omega^2]\sin\beta_0 \\ & + \left[ (\ddot{u}+\ddot{v},z\frac{D_L}{2}) + 2(\dot{v}+\dot{v},z\frac{D_L}{2})\Omega - (u+u,z\frac{D_L}{2})\Omega^2 \right] \cos\beta_0 \end{aligned} \right]_{Z=\ell} \sin(\Omega t + \vartheta_j) \quad (14) $$

The centrifugal stress in the blade at any typical span location 's' from the root of the airfoil is given by,

$$ (\sigma_{cf})_s = (\rho_b \Omega^2 / 2) [L^2 + 2rL - 2rs - s^2] \quad (15) $$

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=(\rho_b cL)$ is expressed as,

$$ F_{cf} = M_b (r+L/2) \Omega^2 \quad (16) $$

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 ($\beta_L - \beta_0$) in the blade, the tendency of which is to untwist the free tip-end of the cantilevered configuration, and is determined as,

$$ T = F_{cf} (C/2) \sin^2 (\beta_L - \beta_0) \quad (17) $$

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,

$$ \sum_{j=1}^{N_b} (F_X)_j = - \sum_{j=1}^{N_b} [(F_{\zeta} - M_b \ddot{\zeta}_0)_s = 0 \cos(\vartheta_j + \Omega t) + (F_{\eta} - M_b \ddot{\eta}_0)_s = 0 \cos\beta_0 \sin(\vartheta_j + \Omega t)]_j \quad (18) $$

and,

$$ \sum_{j=1}^{N_b} (F_Y)_j = - \sum_{j=1}^{N_b} [(F_{\zeta} - M_b \ddot{\zeta}_0)_s = 0 \sin(\vartheta_j + \Omega t) - (F_{\eta} - M_b \ddot{\eta}_0)_s = 0 \cos\beta_0 \cos(\vartheta_j + \Omega t)]_j \quad (19) $$

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,

$$ \sum_{j=1}^{N_b} (F_X)_j = - \sum_{j=1}^{N_b} [(F_{\zeta})_s = 0 \cos(\vartheta_j + \Omega t) + (F_{\eta})_s = 0 \cos\beta_0 \sin(\vartheta_j + \Omega t)]_j \quad (20) $$

and,

$$ \sum_{j=1}^{N_b} (F_Y)_j = - \sum_{j=1}^{N_b} [(F_{\zeta})_s = 0 \sin(\vartheta_j + \Omega t) - (F_{\eta})_s = 0 \cos\beta_0 \cos(\vartheta_j + \Omega t)]_j \quad (21) $$

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 asymmetry 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 $\theta_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 non-symmetric 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 $\theta_0$ in the global frame of reference are:

$$ (BM_U)_X = M_b (r+L/2) \Omega^2 (C/2) \sin^2 (\beta_L - \beta_0) \sin(\theta_0 + \Omega t) \quad (22) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Schematic front view of the bladed-fan-disk with asymmetric rotor due to blade-loss during a typical fan-blade-out (FBO) event.</figcaption></figure>

$$(BM_U)_Y = M_b(r+L/2)\Omega^2(C/2)\sin^2(\beta_L - \beta_0)\cos(\theta_0 + \Omega t) \quad (23)$$

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 $\Omega$ 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,

$$[\zeta_0(t)]_j = u(\ell)\cos(\Omega t + \theta_j) + v(\ell)\sin(\Omega t + \theta_j) + u'(\ell)(D_L/2)\cos(\Omega t + \theta_j) + v'(\ell)(D_L/2)\sin(\Omega t + \theta_j) \quad (24)$$

and,

$$[\eta_0(t)]_j = \cos\beta_0[u(\ell)\sin(\Omega t + \theta_j) - v(\ell)\cos(\Omega t + \theta_j)] + u'(\ell)[\sin\beta_0\cos(\Omega t + \theta_j) + (D_L/2)\cos\beta_0\sin(\Omega t + \theta_j)] + v'(\ell)[\sin\beta_0\sin(\Omega t + \theta_j) - (D_L/2)\cos\beta_0\cos(\Omega t + \theta_j)] \quad (25)$$

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:

$$u(Z,t) = \sum_{i=0}^{\infty} U_i(t)w_i(Z) = U_0(t) + \sum_{i=1}^{\infty} U_i(t)[1 - \cos\phi_i Z] \quad (26)$$

$$v(Z,t) = \sum_{i=0}^{\infty} V_i(t)w_i(Z) = V_0(t) + \sum_{i=1}^{\infty} V_i(t)[1 - \cos\phi_i Z] \quad (27)$$

$$\eta(s,c,t) = \sum_{j=0}^{\infty} \left[ \sum_{i=0}^{\infty} Y_i(t)W_i(s) \right] y_j(c) = Y_0(t) + \sum_{j=1}^{\infty} \left[ \sum_{i=1}^{\infty} Y_i(t)W_i(s) \right] y_j \cos\alpha_j \left(c + \frac{C}{2}\right) \quad (28)$$

$$\zeta(s,t) = \sum_{i=0}^{\infty} Z_i(t)X_i(s) = Z_0(t) + \sum_{i=1}^{\infty} Z_i(t)X_i(s) = Z_0(t) + \sum_{i=1}^{\infty} Z_i(t)[\sin\phi_i s] \quad (29)$$

We define $\phi_i = ((2i-1)\pi)/2\ell$, $\phi_i = (((2i-1)\pi)/2L)\Psi_2$, and, $\alpha_j = (((j-1)\pi)/C)\Psi_1$ with nondimensional terms for the geometrical parameters of the airfoil as $\Psi_1$, $\Psi_2$ (see Appendix B). Then using Eq. (29) and taking its first derivative with respect to s, we also obtain,

$$\zeta_{s}(s,t)|_{s=0} = \sum_{i=1}^{\infty} Z_{i}(t)\phi_{i}\cos\phi_{i}s|_{s=0} = \sum_{i=1}^{\infty} Z_{i}(t)\phi_{i} \quad (30)$$

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**, **C** and **K** matrices. In a gyroscopic system, the velocity-dependent matrix **C** is usually broken into two parts called **C**<sub>D</sub> for the damping effect and **C**<sub>G</sub> 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,

$$ M\ddot{\mathbf{f}}(t) + (\mathbf{C}_D + \mathbf{C}_G)\dot{\mathbf{f}}(t) + (\mathbf{K}_\sigma + \mathbf{K}_F + \mathbf{K}_\Omega)\mathbf{f}(t) = \mathbf{p}(t) \quad (31) $$

Similarly, components of stiffness matrix such as $K_\sigma$ contains an elastic modulus-dependent stiffness, whereas $K_F$ and $K_\Omega$ 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 $\mathbf{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 $\eta(s,c,t)$ in the lateral direction, which results into rigid-body translational inertia mass effect at the airfoil clamped root being as:

$$ M_b \ddot{\eta}_0(t) = M_b \cos\beta_0[\sin(\beta_j + \Omega t) \ddot{U}_0 - \cos(\beta_j + \Omega t) \dot{V}_0] \quad (32) $$

Similarly in the longitudinal direction of the blade due to deformation function $\zeta(s,t)$, the respective rigid-body translational inertia mass effect at the clamped-end becomes:

$$ M_b \ddot{\zeta}_0(t) = M_b[\cos(\beta_j + \Omega t) \ddot{U}_0 + \sin(\beta_j + \Omega t) \dot{V}_0] \quad (33) $$

in addition, the flexural bending of the blades introduces additional inertia terms as,

(a) Blade rotary inertia due to tangential motion:

$$ \rho_b \mathrm{Ch} \int_0^L (r+s)^2 \cos\beta W'_p \mathrm{d}s \quad (34) $$

(b) Blade gyroscopic inertia due to axial motion:

$$ \rho_b \mathrm{Ch} \int_0^L (r+s)^2 \sin\beta W'_p \mathrm{d}s \quad (35) $$

Thus, the global mass and stiffness matrices $\mathbf{M}$ and $\mathbf{K}$ for the current flexible bladed-rotor system can be expressed in terms of multiple sub-matrices. The global stiffness matrix $\mathbf{K}$ is made up of smaller sub-matrices such as, $\mathbf{K}_{(U,U)}$, $\mathbf{K}_{(U,V)}$, $\mathbf{K}_{(V,V)}$, $\mathbf{K}_{(\eta,U_0)}$, $\mathbf{K}_{(\eta,V_0)}$ and $\mathbf{K}_{(\eta,\eta)}$ etc. Each term of the sub-matrices such as $\mathbf{K}_{(\eta,U_0)}$, $\mathbf{K}_{(\eta,V_0)}$ and $\mathbf{K}_{(\eta,\eta)}$ etc. are derived individually and is shown in Appendix B. Using these sub-matrices and the computed values of generalized coordinates like $\mathbf{u}$, $\mathbf{v}$ and $\mathbf{y}^{(j)}$, $\mathbf{z}^{(j)}$ etc. for any typical blade, the force-displacement relationship in the pseudo-static condition can be written as,

$$ \mathbf{K}_{(\eta,U_0)}^{(j)} U_0 + \mathbf{K}_{(\eta,V_0)}^{(j)} V_0 + \mathbf{K}_{(\eta,U)}^{(j)} \mathbf{u} + \mathbf{K}_{(\eta,V)}^{(j)} \mathbf{v} + \mathbf{K}_{(\eta,\eta)}^{(j)} \mathbf{y}^{(j)} + \mathbf{K}_{(\eta,\zeta)}^{(j)} \mathbf{z}^{(j)} = \mathbf{p}_{\eta}^{(j)} \quad (36) $$

$$ \mathbf{K}_{(\zeta,U_0)}^{(j)} U_0 + \mathbf{K}_{(\zeta,V_0)}^{(j)} V_0 + \mathbf{K}_{(\zeta,U)}^{(j)} \mathbf{u} + \mathbf{K}_{(\zeta,V)}^{(j)} \mathbf{v} + \mathbf{K}_{(\zeta,\eta)}^{(j)} \mathbf{y}^{(j)} + \mathbf{K}_{(\zeta,\zeta)}^{(j)} \mathbf{z}^{(j)} = \mathbf{p}_{\zeta}^{(j)} \quad (37) $$

In a steady-state condition of a typical rotating blade no. 'j', the external force term $\mathbf{p}_\eta^{(j)}$ is due to centrifugal forces and $\mathbf{p}_\zeta^{(j)}$ is caused by the blade tip-rub. Using the superscript (j) for the fan blade sequence number (j=1...N), the corresponding rotodynamic equations of motion in the matrix form for the flexible bladed-rotor system under consideration can be

<figure><img src="image_5.png" /><figcaption>Fig. 4. Components of friction forces at the blade-tip during the fan blade rub with the outer casing in the global frame of reference.</figcaption></figure>
written as,

$$ \begin{bmatrix} \mathbf{M}_{S2(m+1) \times 2(m+1)} & 0 & 0 & 0 & 0 \\ 0 & \mathbf{M}_{B}^{(j)} & 0 & 0 & 0 \\ 0 & 0 & \bullet & 0 & 0 \\ 0 & 0 & 0 & \mathbf{M}_{B}^{(N-1)} & 0 \\ 0 & 0 & 0 & 0 & \mathbf{M}_{B}^{(N)} \end{bmatrix} \begin{bmatrix} \ddot{\mathbf{\delta}}_{S2(m+1)} \\ \ddot{\mathbf{\delta}}_{B}^{(j)}(n1+n2) \\ \bullet \\ \ddot{\mathbf{\delta}}_{B}^{(N-1)}(n1+n2) \\ \ddot{\mathbf{\delta}}_{B}^{(N)}(n1+n2) \end{bmatrix} + \tag{38} $$

$$ \begin{bmatrix} \mathbf{K}_{S2(m+1) \times 2(m+1)} & \mathbf{K}_{SB}^{(j)}(2(m+1) \times X(n1+n2)) & \bullet & \mathbf{K}_{BS}^{(N-1)} & \mathbf{K}_{SB}^{(N)} \\ \mathbf{K}_{BS}^{(j)}(n1+n2) \times X(n1+n2) & \mathbf{K}_{B}^{(j)}(n1+n2) & 0 & 0 & 0 \\ \bullet & 0 & \bullet & 0 & 0 \\ \mathbf{K}_{BS}^{(N-1)} & 0 & 0 & \mathbf{K}_{B}^{(N-1)} & 0 \\ \mathbf{K}_{BS}^{(N)} & 0 & 0 & 0 & \mathbf{K}_{B}^{(N)} \end{bmatrix} \begin{bmatrix} \frac{\delta_{S2(m+1)}}{\delta_{B}^{(j)}(n1+n2)} \\ \bullet \\ \frac{\delta_{B}^{(N-1)}(n1+n2)}{\delta_{B}^{(N)}(n1+n2)} \end{bmatrix} = \begin{bmatrix} \frac{\mathbf{p}_{S2(m+1)}}{\mathbf{p}_{B}^{(j)}(n1+n2)} \\ \bullet \\ \frac{\mathbf{p}_{B}^{(N-1)}(n1+n2)}{\mathbf{p}_{B}^{(N)}(n1+n2)} \end{bmatrix} $$

In the above global matrix, the sign ‘●’ is used to indicate the generic form of fully-populated sub-matrices with nonzero terms such as, $\mathbf{K}_{BS}^{(j)}$, $\mathbf{K}_{SB}^{(j)}$, and $\mathbf{K}_{B}^{(j)}$ 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-matricess such as $\mathbf{K}_{BS}^{(j)}$ and $\mathbf{K}_{SB}^{(j)T}$, due to cross-coupling effect between rotor-shaft and the blades, are actually time-dependent with sinusoidal functions containing terms like $\sin(\vartheta_j + \Omega t)$ and $\cos(\vartheta_j + \Omega t)$. It should be obvious that the above sub-matrices pertaining to individual fan blades, are repeated along the diagonal of the global stiffness matrix $\mathbf{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 $\mathbf{M}$ and the velocity-dependent coefficient matrix $\mathbf{C}$. Similarly in a nonsymmetric bladed rotor, we have time-dependent terms appearing in the mass sub-matrix $\mathbf{M}_s$ with parameters like $\sin2(\vartheta_j + \Omega t)$ and $\cos2(\vartheta_j + \Omega t)$.

## 2.2. Displacement and external force vectors for the rotor-shaft due to blade-loss

In Eq. (40), the respective displacement vector terms for the rotor-shaft and the individual blades such as $\delta_s$ and $\delta_B^{(j)}$ in the expanded form are,

$$ \delta_{S2(m+1)}^{(j)} = \begin{Bmatrix} U_0 \\ V_0 \\ \mathbf{u}_m \\ \mathbf{v}_m \end{Bmatrix} \quad \text{and} \quad \delta_{B(n1+n2)}^{(j)} = \begin{Bmatrix} \mathbf{y}_{n1}^{(j)} \\ \mathbf{z}_{n2}^{(j)} \end{Bmatrix} \tag{39} $$

In this section, the terms of the generalized displacement vector $\mathbf{f}(t)$ and the external force vector $\mathbf{p}(t)$ are shown in the expanded form as,

$$ \mathbf{f}(t) = \begin{Bmatrix} U_{\theta}(t) \\ V_{\theta}(t) \\ \mathbf{u}(t)_{m} \\ \mathbf{v}(t)_{m} \\ \bullet \\ \bullet \\ \bullet \\ \bullet \\ \mathbf{y}^{(j)}(t)_{n1} \\ \mathbf{z}^{(j)}(t)_{n2} \\ \bullet \\ \bullet \\ \mathbf{y}^{(N)}(t)_{n1} \\ \mathbf{z}^{(N)}(t)_{n2} \end{Bmatrix} \quad \text{and,} \quad \mathbf{p}(t) = \begin{Bmatrix} F_{X\theta}(t) \\ F_{Y\theta}(t) \\ \mathbf{p}_{X}(t)_{m} \\ \mathbf{p}_{Y}(t)_{m} \\ \bullet \\ \bullet \\ \bullet \\ \bullet \\ \mathbf{p}_{\eta}^{(j)}(t)_{n1} \\ \mathbf{p}_{\zeta}^{(j)}(t)_{n2} \\ \bullet \\ \bullet \\ \mathbf{p}_{\eta}^{(N)}(t)_{n1} \\ \mathbf{p}_{\zeta}^{(N)}(t)_{n2} \end{Bmatrix} \tag{40} $$

The external forces $\mathbf{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 $\mathbf{M}$ are no longer constant rather they oscillate with twice the frequency of the spin. Due to large imbalance, the tips of the fan
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,

$$P(t) = M_b \sum_{j=1}^{N-1} A_{\text{axial}}^{(j)} \quad (48)$$

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 $(1 \le q \le m)$ 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,

$$P(t) = M_b r \sum_{q=1}^{m} w'_{q}(t) \left[ \begin{aligned} & [\ddot{u}_q + 2\dot{v}_q\Omega - u_q\Omega^2]\cos(\theta_0 + \pi + \Omega t) \\ & + [\ddot{v}_q - 2\dot{u}_q\Omega - v_q\Omega^2]\sin(\theta_0 + \pi + \Omega t) \end{aligned} \right] \quad (49)$$

In the lateral force relationship on the blade-tip shown in Eq. (43), the sign of the $\eta$-component of the rub forces are established by keeping it opposite of the sign of $\dot{\eta}(L,c)$, which in the inertial frame of reference is determined as follows:

$$\dot{\eta}(L,c) = \dot{Y}_0(t) + \sum_{l=1}^{n_1} \dot{Y}_l(t) W_l(L) y_p(c) \quad (50)$$

We also consider that the blade-tip rigid body velocity in the tangential direction due to spinning of rotor is $-\Omega R$, whereas the blade clamped-root tangential velocity due to shaft oscillation is determined by,

$$\dot{Y}_0(t) = \left[ \begin{aligned} & \left[ \dot{U}_0 + \sum_{q=1}^{m} \dot{\mathbf{u}}_q(t) \left( w_q(t) + \frac{D_L}{2} w'_q(t) \right) \right] \sin(\theta_j + \Omega t) \\ & - \left[ \dot{V}_0 + \sum_{q=1}^{m} \dot{\mathbf{v}}_q(t) \left( w_q(t) + \frac{D_L}{2} w'_q(t) \right) \right] \cos(\theta_j + \Omega t) \end{aligned} \right] \cos\beta_L \quad (51)$$

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,

$$T(t) = T_0 + T_t \cos \Omega t \quad (52)$$

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)_{\text{New}}=(\Omega_{\text{New}}/\Omega_{\text{Old}})^2(T_0)_{\text{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 $\mu$, the depth of the radial incursion $\epsilon_j(t)$ of the blade-tip inside the abradable surface and the radial stiffness $K_{\text{Rub}}$ given as,

$$T_t = \mu(r+L)K_{\text{Rub}}\epsilon_j(t) \quad (53)$$

The general solution of Eq. (31) can be assumed in the form of a Fourier series with the generalized coordinate $\mathbf{f}(t)$ expressed as,

$$\mathbf{f}(t) = e^{\lambda t} \left[ \frac{\mathbf{b}_0}{2} + \sum_{q=1}^{\infty} \left[ \mathbf{a}_q \sin(q\Omega t) + \mathbf{b}_q \cos(q\Omega t) \right] \right] \quad (54)$$

Then the corresponding velocity and acceleration column vectors $\dot{\mathbf{f}}(t)$ and $\ddot{\mathbf{f}}(t)$ are written as,

$$\dot{\mathbf{f}}(t) = e^{\lambda t} \left[ \frac{\lambda \mathbf{b}_0}{2} + \sum_{q=1}^{\infty} \left[ (\lambda \mathbf{a}_q - q\Omega) \mathbf{b}_q \sin(q\Omega t) + (\lambda \mathbf{b}_q + q\Omega) \mathbf{a}_q \cos(q\Omega t) \right] \right] \quad (55)$$

and,

$$\ddot{\mathbf{f}}(t) = e^{\lambda t} \left[ \frac{\lambda^2 \mathbf{b}_0}{2} + \sum_{q=1}^{\infty} \left[ \begin{aligned} & \left[ (\lambda^2 - (q\Omega)^2) \mathbf{a}_q - 2\lambda q\Omega \mathbf{b}_q \right] \sin(q\Omega t) \\ & + \left[ (\lambda^2 - (q\Omega)^2) \mathbf{b}_q + 2\lambda q\Omega \mathbf{a}_q \right] \cos(q\Omega t) \end{aligned} \right] \right] \quad (56)$$

where, the Fourier coefficient terms such as $\mathbf{b}_0, \mathbf{a}_q, \mathbf{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 $\dot{\phi}$ of the shaft at any axial location can also be determined as,

$$ \dot{\phi} = \frac{d}{dt} \left[ \tan^{-1} \frac{u}{v} \right] \text{ radians/s} \qquad (57) $$

Thus, the shaft whirl frequency in Hz at the axial location of the fan disk is simplified as

$$ \text{Fan disk whirl frequency} = \frac{1}{2\pi[1 + (u/v)^2]} \left[ \frac{\dot{u}}{v} - \frac{u}{v^2} \dot{v} \right]_{Z=\ell} \text{ Hz} \qquad (58) $$

## 2.3. 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 $\delta_s(z_i)$ exceeds the pre-specified value of the corresponding radial gap $\gamma_i$ of the bearing in that location, such that

$$ \delta_s(z_i) = \left[ \sqrt{u(z_i, t)^2 + v(z_i, t)^2} \right]_{\ell} > \gamma_i \qquad (59) $$

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,

$$ (F_X)_i^b = [-S_{xx}^{b(i)}(u-\gamma_i \cos \varphi) - S_{xy}^{b(i)}(v-\gamma_i \sin \varphi) - B_{xx}^{b(i)}u_{,ZZ}]_{Z=z_i} \qquad (60) $$

$$ (F_X)_i^b = [-S_{yx}^{b(i)}(u-\gamma_i \cos \varphi) - S_{yy}^{b(i)}(v-\gamma_i \sin \varphi) - B_{yy}^{b(i)}v_{,ZZ}]_{Z=z_i} \qquad (61) $$

where, the shaft radial deformation angle for a particular bearing located at the axial coordinate of $Z=z_i$ is determined as,

$$ \varphi = \tan^{-1} \frac{u(z_i, t)}{v(z_i, t)} \text{ radians} \qquad (62) $$

# 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,

$$ K_{\text{Rub}} = A_b (E/h)_{\text{Abradable}} \qquad (63) $$

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 $\mu=0.01$.

## 3.1. 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 ... 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.
<table><caption>**Table 1**
Details of the parameters considered for the sample bladed-disk rotor.</caption><tbody><tr><td>M<sub>b</sub></td><td>Mass of one blade (M<sub>b</sub> = ρ<sub>b</sub>ChL)</td><td>10 kg</td></tr><tr><td>(M<sub>D</sub>)<sub>0</sub></td><td>Mass of the overhung fan-disk</td><td>150 kg</td></tr><tr><td></td><td>Mass of additional rotor components (spinner cone, booster spool)</td><td>280 kg</td></tr><tr><td>D<sub>L</sub></td><td>Axial length of the fan-disk (rigid)</td><td>25 cm</td></tr><tr><td>R</td><td>Blade tip radius</td><td>150 cm</td></tr><tr><td>r</td><td>Blade root radius = disk Outer radius</td><td>50 cm</td></tr><tr><td>L</td><td>Blade span length</td><td>100 cm</td></tr><tr><td>A<sub>b</sub></td><td>Blade cross-sectional area</td><td>5000 mm<sup>2</sup></td></tr><tr><td>I<sub>b</sub></td><td>Blade cross-sectional area moment of inertia</td><td>0.5E6 mm<sup>4</sup></td></tr><tr><td>K<sub>Rub</sub></td><td>Case filler material (abradable) stiffness</td><td>10 M N/m</td></tr><tr><td>ℓ</td><td>Length of the shaft</td><td>350 cm</td></tr><tr><td>r<sub>o</sub></td><td>Outside radius of the rotor shaft</td><td>7.0 cm</td></tr><tr><td>r<sub>i</sub></td><td>Inside radius of the rotor shaft</td><td>5.5 cm</td></tr><tr><td>ρ<sub>s</sub></td><td>Mass density of the shaft material</td><td>7.833 g/cc</td></tr><tr><td>β<sub>L</sub></td><td>Blade stagger angle at the tip</td><td>-60°</td></tr><tr><td>Ω</td><td>Initial spin velocity of rotor (radian/s)</td><td>314.1 rad/s</td></tr><tr><td>m<sub>r</sub></td><td>Mass moment of imbalance about engine axis</td><td>1E6 g cm</td></tr><tr><td>θ<sub>0</sub></td><td>Circumferential location of the unbalance at time 't'=0 (constant)</td><td>180°</td></tr><tr><td>N</td><td>Original number of blades on the fan-rotor</td><td>20</td></tr><tr><td>N<sub>c</sub></td><td>Number of discrete bearings supporting the shaft</td><td>3</td></tr><tr><td>Z<sub>i</sub></td><td>Bearing support locations</td><td>0.3, 266.7 cm, 325 cm</td></tr><tr><td>γ</td><td>Blade radial tip clearance (nonrotating)</td><td>2 cm</td></tr><tr><td>μ</td><td>Coefficient of friction between the blade tip and abradable (dynamic)</td><td>0.01</td></tr><tr><td>E<sub>s</sub></td><td>Elastic Young's modulus of the rotor shaft material</td><td>200 G Pa</td></tr><tr><td>E<sub>b</sub></td><td>Elastic Young's modulus of the blade material</td><td>117 G Pa</td></tr><tr><td>C<sub>0</sub></td><td>Blade chord at the root</td><td>25 cm</td></tr><tr><td>C<sub>L</sub></td><td>Blade chord at the tip</td><td>50 cm</td></tr><tr><td>C</td><td>Average chord of the blade airfoil</td><td>37.5 cm</td></tr><tr><td>Rotational speed (rev/min)</td><td></td><td>3000 rev/min</td></tr><tr><td>Roller Bearing-1: at Z=z<sub>1</sub>=0.3 cm</td><td>S<sub>xx</sub><sup>b(1)</sup> = S<sub>yy</sub><sup>b(1)</sup> = 2000.0 M N/m</td><td></td></tr><tr><td>Roller Bearing-2: at Z=z<sub>2</sub>=265.0 cm</td><td>S<sub>xx</sub><sup>b(2)</sup> = S<sub>yy</sub><sup>b(2)</sup> = 715.0 M N/m</td><td></td></tr><tr><td>Ball Bearing-3: at Z=z<sub>3</sub>=325.0 cm</td><td>S<sub>xx</sub><sup>b(3)</sup> = S<sub>yy</sub><sup>b(3)</sup> = 565.0 M N/m and B<sub>xx</sub><sup>b(3)</sup> = B<sub>yy</sub><sup>b(3)</sup> = 7.500E6 N m/radian</td><td></td></tr><tr><td>Bearing-3 Radial load capability</td><td>1 M N</td><td></td></tr></tbody></table>

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 Hz bifurcates in two vibrational modes with the lowest rotor modes showing at 58.35 Hz and the next one slightly higher at 59.97 Hz. Similarly, the other mode previously showing up at 70.26 Hz shifts along 70.26 Hz and 70.28 Hz, respectively. In the rotating condition with spinning speed as 1000 rev/min and 3000 rev/min, the corresponding natural frequencies of this rotor are computed as follows:

(a) 61.31 Hz and 70.30 Hz at 1000 rev/min (rotating condition)
62.85 Hz and 70.35 Hz at 1000 rev/min (rotating condition)

(b) 63.23 Hz and 70.20 Hz at 3000 rev/min (rotating condition)
64.44 Hz and 70.21 Hz at 3000 rev/min (rotating condition).

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.
<figure><img src="image_6.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Fan blade frequencies and respective mode shapes in stationary and running conditions: (a) 1F-mode (45.89 Hz stationary blade, 51.65 Hz at 1000 rev/min, 84.08 Hz at 3000 rev/min), (b) 1T-mode (182.78 Hz stationary blade, 183.44 Hz at 1000 rev/min, 198.89 Hz at 3000 rev/min) and (c) 2T-mode (552.50 Hz stationary blade, 550.82 Hz at 1000 rev/min, 549.86 Hz at 3000 rev/min).</figcaption></figure>

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 drops to 39.67 Hz with only Bearings-1 and -2 remaining active as a result of the fusing of Bearing-3 soon after encountering the large imbalance caused by the fan-blade-loss. The shaft mode frequencies of 39.67 Hz, sometime
<figure><img src="image_9.png" /><figcaption>Fig. 6. Symmetric rotor natural modes with rigid blades (no blade participation).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Symmetric rotor natural modes with flexible blades (blades deforming in its own flexural vibrational mode). (a) 58.45 Hz, And 70.26 Hz at zero speed (stationary condition) and (b) 62.73 Hz and 70.20 Hz at 3000 rev/min (rotating condition).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Critical speed determination by the forced harmonic response analysis (u,v) of the rotor at the bearing support location for a nominal rotating unbalance as the speed drops.</figcaption></figure>
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 blade-tips 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.

## 3.2. 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(t)$ and $v(t)$ 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 ($\eta_i, \zeta_i$), 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 $\pm 2$ 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
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Analytically predicted rotor orbits of an asymmetric rotor-shaft at the fan disk-center subsequent to a sudden blade-loss at 3000 rev/min.
(a) Rigid blade formulation and (b) deformable blade formulation.</figcaption></figure>

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 rigid-
blade 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.
<figure><img src="image_14.png" /><figcaption>Fig. 10. Radial excursion of the rotor shaft at the fan disk-center after blade-loss at 3000 rev/min (deformable versus rigid blade solution).</figcaption></figure>

The current pre-twisted blade model automatically takes into account of many different types of complex mode shapes as well (see Fig. 5).

## 4. 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:

(a) Horizontal and vertical movement of the fan disk.

(b) Blade tip-rub load time-history.

(c) Rotor response frequency.

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.
<figure><img src="image_15.png" /><figcaption>Fig. 11. Time-domain horizontal response of a nonsymmetric shaft at the fan disk-center during windmilling at 1000 rev/min after a blade-loss.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Zoomed-in view of a typical 1 s window of the horizontal response of the nonsymmetric shaft windmilling at 1000 rev/min after a blade-loss.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Time-domain vertical response of a nonsymmetric shaft at the fan disk-center during windmilling at 1000 rev/min after a blade-loss.</figcaption></figure>

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 $Ω/2$, and the system responds at its sum-frequency of 25.4 Hz also.
<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Zoomed-in view of a typical 1 s window of the vertical response of the nonsymmetric shaft windmilling at 1000 rev/min after a blade-loss.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Time-domain response of the radial rub load on the rotor of the asymmetric shaft at the fan disk-center during windmilling after a blade-loss.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16.** Frequency content of vibratory response of a windmilling nonsymmetric fan rotor windmilling at 1000 rev/min (16.66 Hz) at the bearing support.</figcaption></figure>

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.
<figure><img src="image_21.png" /><figcaption>Fig. 17. Whirling frequencies of a windmilling nonsymmetric fan rotor at 1000 rev/min (16.66 Hz) with a blade missing.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Comparison of measured airfoil (L=4 cm, $L/C=10517$, $h=0.125$ cm, $\beta=-45^\circ$) root strain gage data (...) versus results from the analytical blade-tip rub model (—) during repeated radial incursion of 0.1 mm at the blade-tip with 72° circumferential rub-zone for the first 4 rub-events.</figcaption></figure>

This rotor depending upon the blade tip clearance $\gamma$ 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-dependent circumferential location of a typical blade identified as 'Blade-j' which might be rubbing at any time 't' is expressed as, $(\theta_j+\Omega t)$. 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].

## 5. 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.

## 6. 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 chord-locations, 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.

## Acknowledgments

The author would like to acknowledge Mr. Sujit Ojha and Mr. Illya Arcos for their assistance in preparing the manuscript of this paper. The author would also thank his manager Mr. Mohamed Ali for reviewing the draft version of the paper and his constructive feed-back as well as for supporting its presentation.

## Appendix A. Rotordynamic equations of motion of a nonsymmetric rotor

### (A1) Gyroscopic effect of individual fan blades from energy formulation

The angular instantaneous velocity of any typical j-th fan blade at its mass-center due to oscillatory motion of the shaft-disk assembly in the global frame of reference is expressed as, $\omega = \dot{\theta}_X \hat{i} + \dot{\theta}_Y \hat{j}$. The instantaneous position vector of the mass-center of that blade with respect to the origin at the engine axis is located at,

$$ \mathbf{r}_{\text{CG}} = (r+L/2)[\cos(\vartheta_j + \Omega t)\hat{i} + \sin(\vartheta_j + \Omega t)\hat{j}] \quad (\text{A1}) $$

The axial movement in the global Z-direction of any typical jth blade-CG due to rigid-body rotation of the disk is:

$$ \theta_X(r+L/2)\sin(\vartheta_j + \Omega t) - \theta_Y(r+L/2)\cos(\vartheta_j + \Omega t) \quad (\text{A2}) $$

which yields the axial velocity vector of the blade mass center (C.G.) in the global frame of reference as,

$$ = (r+L/2)[[\Omega u_z - \dot{v}_z]\sin(\vartheta_j + \Omega t) - [\Omega v_z + \dot{u}_z]\cos(\vartheta_j + \Omega t)]\hat{k} \quad (\text{A3}) $$

Thus, the axial kinetic energy $\Gamma$ of the blade due to the rigid-body oscillatory motion of the fan-disk will be,

$$ (\Gamma_{\text{Axial-KE}})_j = (1/2) M_b (r+L/2)^2 \left[ \begin{array}{l} [\Omega u_z - \dot{v}_z]\sin(\vartheta_j + \Omega t) \\ -[\Omega v_z + \dot{u}_z]\cos(\vartheta_j + \Omega t) \end{array} \right]^2 \quad (\text{A4}) $$

At this point, using the kinetic energy term $\Gamma$, and variational calculus in terms of generalized coordinates $u_z(Z,t)$ and $v_z(Z,t)$ as the displacement variables, we obtain the contribution of blade inertia parameters in forming the general mass and circulatory effects in the equations of motion:

$$ -\frac{d}{dZ} \left[ \frac{d}{dt} \left( \frac{\partial \Gamma}{\partial \dot{u}_z} \right) - \frac{\partial \Gamma}{\partial u_z} \right] + \frac{d}{dt} \left( \frac{\partial \Gamma}{\partial \dot{v}_z} \right) - \frac{\partial \Gamma}{\partial v} = 0 \quad (\text{A5}) $$

and,

$$ -\frac{d}{dZ} \left[ \frac{d}{dt} \left( \frac{\partial \Gamma}{\partial \dot{v}_z} \right) - \frac{\partial \Gamma}{\partial v_z} \right] + \frac{d}{dt} \left( \frac{\partial \Gamma}{\partial \dot{u}_z} \right) - \frac{\partial \Gamma}{\partial u_z} = 0 \quad (\text{A6}) $$
where the time-dependent gyroscopic and inertia terms contributed by each individual blade ($j=1...N_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:

$$ = \sum_{j=1}^{N_b} \frac{M_b(r+L/2)^2}{2} \left[ \begin{aligned} & (\ddot{v}_Z - 2\Omega\dot{u}_Z - \Omega^2 v_Z) - (\ddot{v}_Z - 2\Omega\dot{u}_Z - \Omega^2 v_Z)\cos2(\vartheta_j + \Omega t) \\ & + (\ddot{u}_Z + 2\Omega\dot{v}_Z - \Omega^2 u_Z)\sin2(\vartheta_j + \Omega t) \end{aligned} \right]_{Z=\ell} \quad (A7) $$

and,

(b) Bladed-disk dynamic moment at the shaft center about the global Y-axis:

$$ = \sum_{j=1}^{N_b} \frac{M_b(r+L/2)^2}{2} \left[ \begin{aligned} & (\ddot{u}_Z + 2\Omega\dot{v}_Z - \Omega^2 u_Z) + (\ddot{u}_Z + 2\Omega\dot{v}_Z - \Omega^2 u_Z)\cos2(\vartheta_j + \Omega t) \\ & + (\ddot{v}_Z - 2\Omega\dot{u}_Z - \Omega^2 v_Z)\sin2(\vartheta_j + \Omega t) \end{aligned} \right]_{Z=\ell} \quad (A8) $$

## Appendix B. Determination of terms in conventional M, C and K matrices

For brevity, we introduce coefficients such as $\alpha_i$, $\phi_i$, $\phi_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,

$$ \left. \begin{aligned} \alpha_i &= \frac{(i-1)\pi}{C} \Psi_1 \\ \varphi_i &= \frac{(2i-1)\pi}{2L} \Psi_2 \\ \phi_i &= \frac{(2i-1)\pi}{2\ell} \end{aligned} \right\} \quad (B1) $$

where nondimensional terms are,

$$ \Psi_1 = \left[ 1 + \left(\frac{h}{L}\right) \left(\frac{C}{h}\right)^2 \frac{(\beta'L)^2}{8} \right] $$

and

$$ \Psi_2 = \left[ 1 + \left(\frac{h}{L}\right)^3 \left(\frac{L}{C}\right)^2 \right] $$

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$):

$$ \sum_{i=1}^{m} w_i(Z) = \sum_{i=1}^{m} (1 - \cos\phi_i Z) \quad \text{with } (0 \le Z \le \ell) \quad (B2) $$

(b) Blade lateral $\eta(s,c)$: $\sum_{j=1}^{n_1} [\sum_{i=1}^{n_1} W_i(s)]_j \cos\chi_j(c + (C/2))$ with $(0 \le s \le L)$ and, $(-C/2 \le c \le C/2)$, where

$$ W_i(s) = 1 + c_2 s^2 + c_3 s^3 - \cos(\phi_i s) \quad (B3) $$

(c) Blade longitudinal $\zeta(s)$:

$$ \sum_{i=1}^{n_2} X_i(s) = \sum_{i=1}^{n_2} \sin(\phi_i s) \quad \text{with } (0 \le s \le L) \quad (B4) $$

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''_i(L) = 0$ and, $W'''_i(L) = 0$ in Eq. (B3), one obtains, $c_2 = -3c_3L - (\phi_i^2)\cos(\phi_i L)/2$, along with $c_3 = [(\phi_i^3)\sin(\phi_i L)]/6$. Thus, the term $W_i(s)$ for the lateral deflection of the blade with all boundary conditions satisfied is,

$$ W_i(s) = \left[ 1 - \frac{\phi_i^2}{2} [L\phi_i\sin(\phi_i L) + \cos(\phi_i L)] s^2 + \frac{(\phi_i^3)\sin(\phi_i L)}{6} s^3 - \cos(\phi_i s) \right] \quad (B5) $$
Hence, the transient transverse deformation term $η=η(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,

$$ \eta(s,c,t) = \sum_{j=1}^{4} \left[ \sum_{i=1}^{n_1} Y_i(t)W_i(s) \right]_j \cos\alpha_j \left(c + \frac{C}{2}\right) \quad (B6) $$

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 $ζ=ζ(s,t)$ of any typical blade is expressed with n2 number of terms of the series as,

$$ \zeta(s,t) = \sum_{k=1}^{n_2} Z_k(t)X_k(s) = \sum_{k=1}^{n_2} Z_k(t)\sin(\phi_k s) \quad (B7) $$

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(θ_j + Ωt)$ and $cos(θ_j + Ωt)$ components. However for specific turbofan blade-loss applications, we can write $Ωt=0$, which enables us to simplify the corresponding vectors in terms of $sin θ_j$ and $cos θ_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 (η-degree of freedom) deformation and

$n_2$ = number of deflection shape function terms used to define blade longitudinal (ζ-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}$

$$ [\mathbf{M}_{(U0,U0)}]_{1,1} = \rho_s \int_0^\ell A_s \mathrm{d}Z + \sum_{k=1}^{N_D} (M_D)_k + \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (B8) $$

$$ [\mathbf{M}_{(U0,U)}]_{1,q} = [\mathbf{M}_{(U,U0)}]_{q,1} = \rho_s \int_0^\ell A_s w_q \mathrm{d}Z + \sum_{k=1}^{N_D} w_q(Z_k)(M_D)_k + [w_q(\ell) + (D_L/2)w'_q(\ell)] \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (B9) $$

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=1...N_b$ times for each blade:

$$ [\mathbf{M}_{(U0,\eta)}^{(j)}]_{1,q} = [\mathbf{M}_{(\eta,U0)}^{(j)}]_{q,1} = -M_b((r+L/2)/r)\cos\beta_0\sin(\theta_j + \Omega t) \quad (B10) $$

$$ [\mathbf{M}_{(U0,ζ)}^{(j)}]_{1,q} = [\mathbf{M}_{(ζ,U0)}^{(j)}]_{q,1} = -M_b((r+L/2)/r)\cos(\theta_j + \Omega t) \quad (B11) $$

$$ [\mathbf{M}_{(V0,V0)}]_{2,2} = \rho_s \int_0^\ell A_s \mathrm{d}Z + \sum_{k=1}^{N_D} (M_D)_k + \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (B12) $$

$$ [\mathbf{M}_{(V0,V)}]_{2,q} = [\mathbf{M}_{(V,V0)}]_{q,2} = \rho_s \int_0^\ell A_s w_q \mathrm{d}Z + \sum_{k=1}^{N_D} w_q(Z_k)(M_D)_k + [w_q(\ell) + (D_L/2)w'_q(\ell)] \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (B13) $$

The following 2 sub-matrices of Eqs. (B14) and (B15) are used in the consistent mass matrix formulation and are repeated $j$=number of fan blades ($N_b$) times:

$$ [\mathbf{M}_{(V0,η)}^{(j)}]_{2,q} = [\mathbf{M}_{(η,V0)}^{(j)}]_{q,2} = M_b((r+L/2)/r)\cos\beta_0\cos(\theta_j + \Omega t) \quad (B14) $$
$$[\mathbf{M}_{(V_0, \zeta)}^{(j)}]_{2,q} = [\mathbf{M}_{(\zeta, V_0)}^{(j)}]_{q,2} = -M_b \left( \frac{r+L/2}{r} \right) \sin(\theta_j + \Omega t) \quad (B15)$$

$$[\mathbf{M}_{(U,U_0)}]_{p,1} = \rho_s \int_0^\ell A_s w_p dZ + \sum_{k=1}^{N_D} w_p(Z_k)(M_D)_k + [w_p(\ell) + (D_L/2)w'_p(\ell)] \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (B16)$$

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 $[\mathbf{M}_{(U,U)}]_{p,q}$ time-dependent resulting into twice of the running speed frequency excitation due to contribution of the $\cos^2(\theta_j + \Omega t)$ term.

$$\begin{aligned} [\mathbf{M}_{(U,U)}]_{p,q} = \rho_s \int_0^\ell A_s w_p w_q dZ + (\rho_s/2) \int_0^\ell A_s (r_o^2 + r_i^2) w'_p w'_q dZ &+ \sum_{k=1}^{N_D} [(M_D)_k w_p(Z_k) w_q(Z_k) + 2(J_D)_k w'_p(Z_k) w'_q(Z_k)]_k \\ &+ \left[ (M_D)_0 + M_b \sum_{j=1}^{N_b} j \right] [w_p(\ell) + (D_L/2) w'_p(\ell)][w_q(\ell) + (D_L/2) w'_q(\ell)] \\ &+ \left[ 2(J_D)_0 + M_b (r+L/2)^2 \sum_{j=1}^{N_b} \cos^2(\theta_j + \Omega t) \right] w'_p(\ell) w'_q(\ell) \end{aligned} \quad (B17)$$

$$[\mathbf{M}_{(U,V)}]_{p,q} = [\mathbf{M}_{(V,U)}]_{q,p} = \left[ M_b (r+L/2)^2 \sum_{j=1}^{N_b} [\sin(\theta_j + \Omega t) \cos(\theta_j + \Omega t)] \right] w'_p(\ell) w'_q(\ell) \quad (B18)$$

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:

$$\begin{bmatrix} \mathbf{M}_{(U,\eta)}^{(j)} \\ \mathbf{M}_{(\eta,U)}^{(j)} \end{bmatrix}_{p,q} = \begin{bmatrix} \mathbf{M}_{(\eta,U)}^{(j)} \\ \mathbf{M}_{(\eta,U)}^{(j)} \end{bmatrix}_{q,p} = 0 = -(M_b/L) \begin{bmatrix} [w_p(\ell) + (D_L/2)w'_p(\ell)][\int_0^L W_q \cos \beta ds] \sin(\theta_j + \Omega t) \\ [w'_p(\ell)][\int_0^L (r+s)^2 W'_q \cos \beta ds] \sin^2(\theta_j + \Omega t) \end{bmatrix} \quad (B19)$$

$$\begin{bmatrix} \mathbf{M}_{(U,\zeta)}^{(j)} \\ \mathbf{M}_{(\zeta,U)}^{(j)} \end{bmatrix}_{p,q} = \begin{bmatrix} \mathbf{M}_{(\zeta,U)}^{(j)} \\ \mathbf{M}_{(\zeta,U)}^{(j)} \end{bmatrix}_{q,p} = 0 = -(M_b/L)[w_p(\ell) + (D_L/2)w'_p(\ell)][\int_0^L X_q ds] \cos(\theta_j + \Omega t) \quad (B20)$$

$$[\mathbf{M}_{(V,V_0)}]_{p,2} = \rho_s \int_0^\ell A_s w_p dZ + \sum_{k=1}^{N_D} w_p(Z_k)(M_D)_k + [w_p(\ell) + (D_L/2)w'_p(\ell)] \left[ (M_D)_0 + \sum_{j=1}^{N_b} (M_b)_j \right] \quad (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):

$$\begin{aligned} [\mathbf{M}_{(V,V)}]_{p,q} = \rho_s \int_0^\ell A_s w_p w_q dZ + (\rho_s/2) \int_0^\ell A_s (r_o^2 + r_i^2) w'_p w'_q dZ &+ \sum_{k=1}^{N_D} [(M_D)_k w_p(Z_k) w_q(Z_k) + 2(J_D)_k w'_p(Z_k) w'_q(Z_k)]_k \\ &+ \left[ (M_D)_0 + M_b \sum_{j=1}^{N_b} j \right] [w_p(\ell) + (D_L/2) w'_p(\ell)][w_q(\ell) + (D_L/2) w'_q(\ell)] \\ &+ \left[ 2(J_D)_0 + M_b (r+L/2)^2 \sum_{j=1}^{N_b} \sin^2(\theta_j + \Omega t) \right] w'_p(\ell) w'_q(\ell) \end{aligned} \quad (B22)$$

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:

$$\begin{bmatrix} \mathbf{M}_{(V,\eta)}^{(j)} \\ \mathbf{M}_{(\eta,V)}^{(j)} \end{bmatrix}_{p,q} = \begin{bmatrix} \mathbf{M}_{(\eta,V)}^{(j)} \\ \mathbf{M}_{(\eta,V)}^{(j)} \end{bmatrix}_{q,p} = 0 = (M_b/L) \begin{bmatrix} [w_p(\ell) + (D_L/2)w'_p(\ell)][\int_0^L W_q \cos \beta ds] \cos(\theta_j + \Omega t) \\ [w'_p(\ell)][\int_0^L (r+s)^2 W'_q \cos \beta ds] \cos^2(\theta_j + \Omega t) \end{bmatrix} \quad (B23)$$

$$\begin{bmatrix} \mathbf{M}_{(V,\zeta)}^{(j)} \\ \mathbf{M}_{(\zeta,V)}^{(j)} \end{bmatrix}_{p,q} = \begin{bmatrix} \mathbf{M}_{(\zeta,V)}^{(j)} \\ \mathbf{M}_{(\zeta,V)}^{(j)} \end{bmatrix}_{q,p} = 0 = -(M_b/L)[w_p(\ell) + (D_L/2)w'_p(\ell)][\int_0^L X_q ds] \sin(\theta_j + \Omega t) \quad (B24)$$

$$\begin{aligned} [\mathbf{M}_{(\eta,\eta)}^{(j)}]_{p,q} = (M_b/L) \int_0^L W_p W_q ds - \frac{(M_b/L)h^2}{12} \int_0^L [W_p W_q'' - (\alpha_k^2) W_p W_q] ds \\ + (M_b/L) (\beta')^2 \left[ 1 + v_b \frac{C^2}{12} \left(\frac{C}{L}\right)^2 \alpha_k^2 \right] \int_0^L s^2 W_p W_q ds \end{aligned} \quad (B25)$$
$$[\mathbf{M}_{(\zeta, U0)}^{(j)}]_{p,1} = -M_b((r+L/2)/r)\cos(\vartheta_j + \Omega t) \quad (B26)$$

$$[\mathbf{M}_{(\zeta, V0)}^{(j)}]_{p,2} = -M_b((r+L/2)/r)\sin(\vartheta_j + \Omega t) \quad (B27)$$

$$[\mathbf{M}_{(\zeta, U)}^{(j)}]_{p,q} = -(M_b/L)[w_q(\ell) + (D_L/2)w'_q(\ell)] \left[ \int_0^L X_p ds \right] \cos(\vartheta_j + \Omega t) \quad (B28)$$

$$[\mathbf{M}_{(\zeta, V)}^{(j)}]_{p,q} = -(M_b/L)[w_q(\ell) + (D_L/2)w'_q(\ell)] \left[ \int_0^L X_p ds \right] \sin(\vartheta_j + \Omega t) \quad (B29)$$

$$[\mathbf{M}_{(\zeta, \zeta)}^{(j)}]_{p,q} = (M_b/L) \int_0^L X_p X_q ds \quad (B30)$$

## (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,

$$C = C_G + C_D \quad (B31)$$

The components of gyroscopic matrix term $C_G$ are always skew-symmetric in nature and contain coefficients like $\pm 2\Omega M$ to account for the spinning rotor with,

$$C_{Gp,q} = -C_{Gq,p} \quad (B32)$$

The components of velocity-dependent coefficient matrix term $C_D$ contain the effect of bearing damping as well as any material damping present in the shaft. Most of the times, these components are symmetric in nature, except when the bearing damping values are different in the horizontal and vertical directions such that its $D_{xy}^{(i)} \neq D_{yx}^{(i)}$. Some of the typical nonzero terms of the $C$ matrix are as follows:

$$[C_{G(U0,\eta)}^{(j)}]_{1,q} = -[C_{G(\eta,U0)}^{(j)}]_{q,1} = 2M_b((r+L/2)/r)\cos\beta_0\Omega\cos(\vartheta_j + \Omega t) \quad (B33)$$

$$[C_{(U,U)}]_{p,q} = 2\psi H \left[ \int_0^\ell w_p w_q'''' dZ \right] + \sum_{i=1}^{N_C} D_{xx}^{(i)} w_p(z_i) w_q(z_i) - 2\Omega M_b (r+L/2)^2 \sum_{j=1}^{N_b} w'_p(\ell) w'_q(\ell) \sin(\vartheta_j + \Omega t) \cos(\vartheta_j + \Omega t) \quad (B34)$$

$$[C_{G(U,V)}]_{p,q} = -2\Omega \int_0^\ell I_d w_p w_p'' dZ + 2\Omega (J_D)_0 w_p'(\ell) w_q'(\ell) + 2\Omega M_b (r+L/2)^2 \sum_{j=1}^{N_b} w_p'(\ell) w_q'(\ell) \cos^2(\vartheta_j + \Omega t) \quad (B35)$$

$$[C_{G(V,U)}]_{p,q} = 2\Omega \int_0^\ell I_d w_p w_q'' dZ + 2\Omega (J_D)_0 w_p'(\ell) w_p'(\ell) - 2\Omega M_b (r+L/2)^2 \sum_{j=1}^{N_b} w_p'(\ell) w_q'(\ell) \sin^2(\vartheta_j + \Omega t) \quad (B36)$$

$$[C_{G(V,\eta)}^{(j)}]_{p,q} = -[C_{G(\eta,V)}^{(j)}]_{q,p} = 2(M_b/L)\Omega w_p'(\ell) \left[ \int_0^L (r+s)^2 W'_q \sin \beta ds \right] \sin(\vartheta_j + \Omega t) \cos(\vartheta_j + \Omega t) \quad (B37)$$

$$[C_{G(\zeta,\eta)}^{(j)}]_{p,q} = 2(M_b/L)\Omega \left[ \int_0^L X_p W_q \cos \beta ds \right] \quad (B38)$$

## (B3) Stiffness terms of displacement-dependent sub-matrices $K_{p,q}$

In the following matrices, the $K_{Rub}$ effects are nonlinear as it depends upon the amount of tip-deformation in the dynamical system. Similarly, for the special case of fan blades being considered as rigid, the contact tip-rub stiffness terms is replaced by another term called $K_{Eq-Rub}$. Here, the term $K_{Eq-Rub}$ is defined as an equivalent contact stiffness term, which accounts for the combined effect of blade flexibility as well as the abradable radial rub stiffness. Its contributions are nonzero only if the fan-shaft-disk radial movement $\delta_s$ is large enough to consume the running blade tip-clearance $\gamma$ with the inner surface of the casing filler material, which may result into intermittent tip-rub-induced vibration of the blade. Additionally, the contributions of frictional force due to $\pm \mu K_{Rub}$ effects are to generate loads in the tangential direction of the rotor which give rise to skew-symmetric terms in the stiffness matrix $K$. The instantaneous positive (+) or negative (-) sign for any new time-step $\ell + \Delta t$ is determined based upon the direction of sliding velocity at the previous time-step $t$. It should be further noted that the stiffness matrix $K$ for a rotor spinning with angular velocity $\Omega$, there will be additional
contribution from stress-softening terms determined as $-\Omega^2 M$. Usually, in majority of the rotordynamic problems this effect of stress-softening is negligible, but in a large jet engine fan-blade rotor this effect can be significant and should be accounted for

$$[\mathbf{K}_{(U0,U0)}]_{1,1} = \sum_{i=1}^{N_c} S_{xx}^{b(i)} + K_{\text{Eq-Rub}} \quad (B39)$$

$$[\mathbf{K}_{(U0,V0)}]_{1,2} = \sum_{i=1}^{N_c} S_{xy}^{b(i)} + \mu K_{\text{Eq-Rub}} \quad (B40)$$

$$[\mathbf{K}_{(U0,U)}]_{1,q} = \sum_{i=1}^{N_c} \left[ \begin{aligned} & S_{xx}^{b(i)}[w_q(z_i) + z_i w'_q(z_i)] + B_{yy}^{b(i)} w''_q(z_i) \\ & - (E_s I) w'''_q(z_i) \end{aligned} \right] - (E_s I) w'''_q(\ell) - P(t) w'_q(\ell) + K_{\text{Eq-Rub}} w_q(\ell) \quad (B41)$$

$$[\mathbf{K}_{(U0,V)}]_{1,q} = \sum_{i=1}^{N_c} S_{xy}^{b(i)}[w_q(z_i) + z_i w'_q(z_i)] + \mu K_{\text{Eq-Rub}} w_q(\ell) \quad (B42)$$

In the following 2 stiffness sub-matrices of Eqs. (B43) and (B44), the contributions due to individual blade deflection-dependent terms such as $X_q(0)$, $W_q(0)$ and its derivatives as $X_q(0)$, $W''_q(0)$ are set to zero for a special case, when the assumed deformation-shape functions arbitrarily satisfy the airfoil clamped-end conditions:

$$[\mathbf{K}_{(U0,\eta)}^{(j)}]_{1,q} = [\mathbf{K}_{(\eta,U0)}^{(j)}]_{q,1} = - \sum_{i=1}^{N_c} [S_{xx}^{b(i)} \sin(\theta_j + \Omega t) - S_{xy}^{b(i)} \cos(\theta_j + \Omega t)] \cos\beta_0 \\ + [[CDW''_{q}(0)] \cos\beta_0 \pm \mu K_{\text{Rub}} W_q(L) \cos\beta_L] \sin(\theta_j + \Omega t) \quad (B43)$$

$$[\mathbf{K}_{(U0,\zeta)}^{(j)}]_{1,q} = [\mathbf{K}_{(\zeta,U0)}^{(j)}]_{q,1} = - \sum_{i=1}^{N_c} [S_{xx}^{b(i)} \cos(\theta_j + \Omega t) + S_{xy}^{b(i)} \sin(\theta_j + \Omega t)] - \left[ \frac{E_b C h}{(1-v_b^2)} X'_q(0) + K_{\text{Rub}} X_q(L) \right] \cos(\theta_j + \Omega t) \quad (B44)$$

$$[\mathbf{K}_{(V0,U0)}]_{2,1} = \sum_{i=1}^{N_c} S_{yx}^{b(i)} - \mu K_{\text{Eq-Rub}} \quad (B45)$$

$$[\mathbf{K}_{(V0,V0)}]_{2,2} = \sum_{i=1}^{N_c} S_{yy}^{b(i)} + K_{\text{Eq-Rub}} \quad (B46)$$

$$[\mathbf{K}_{(V0,U)}]_{2,q} = \sum_{i=1}^{N_c} S_{yx}^{b(i)}[w_q(z_i) + z_i w'_q(z_i)] - \mu K_{\text{Eq-Rub}} w_q(\ell) \quad (B47)$$

$$[\mathbf{K}_{(V0,V)}]_{2,q} = \sum_{i=1}^{N_c} \left[ \begin{aligned} & S_{yy}^{b(i)}[w_q(z_i) + z_i w'_q(z_i)] + B_{xx}^{b(i)} w''_q(z_i) \\ & - (E_s I) w'''_q(z_i) \end{aligned} \right] - (E_s I) w'''_q(\ell) - P(t) w'_q(\ell) + K_{\text{Eq-Rub}} w_q(\ell) \quad (B48)$$

Using the properties of geometry-related deformation-shape functions outlined earlier related to Eqs. (B43), (B44), in the following Eqs. (B49) and (B50), the contributions due to individual blade (j)'s deflection-dependent terms such as $X_q(0)$, $W_q(0)$ and its derivatives as $X_q(0)$, $W''_q(0)$, can again be set to zero:

$$[\mathbf{K}_{(V0,\eta)}^{(j)}]_{2,q} = [\mathbf{K}_{(\eta,V0)}^{(j)}]_{q,2} = - \sum_{i=1}^{N_c} [S_{yx}^{b(i)} \sin(\theta_j + \Omega t) - S_{yy}^{b(i)} \cos(\theta_j + \Omega t)] \cos\beta_0 \\ - [[CDW''_{q}(0)] \cos\beta_0 \pm \mu K_{\text{Rub}} W_q(L) \cos\beta_L] \cos(\theta_j + \Omega t) \quad (B49)$$

$$[\mathbf{K}_{(V0,\zeta)}^{(j)}]_{2,q} = [\mathbf{K}_{(\zeta,V0)}^{(j)}]_{q,2} = - \sum_{i=1}^{N_c} [S_{yx}^{b(i)} \cos(\theta_j + \Omega t) + S_{yy}^{b(i)} \sin(\theta_j + \Omega t)] - \left[ \frac{E_b C h}{(1-v_b^2)} X'_q(0) + K_{\text{Rub}} X_q(L) \right] \sin(\theta_j + \Omega t) \quad (B50)$$

$$[\mathbf{K}_{(U,U0)}]_{p,1} = \sum_{i=1}^{N_c} \left[ \begin{aligned} & S_{xx}^{b(i)}[w_p(z_i) + z_i w'_p(z_i)] + B_{yy}^{b(i)} w''_p(z_i) \\ & - (E_s I) w'''_p(z_i) \end{aligned} \right] - (E_s I) w'''_p(\ell) - P(t) w'_p(\ell) + K_{\text{Eq-Rub}} w_p(\ell) \quad (B51)$$

$$[\mathbf{K}_{(U,V0)}]_{p,2} = \sum_{i=1}^{N_c} S_{xy}^{b(i)}[w_p(z_i) + z_i w'_p(z_i)] + \mu K_{\text{Eq-Rub}} w_p(\ell) \quad (B52)$$

$$[\mathbf{K}_{(U,U)}]_{p,q} = \sum_{i=1}^{N_c} [S_{xx}^{b(i)} w_p(z_i) w_q(z_i) + B_{yy}^{b(i)} w'_p(z_i) w''_q(z_i)] + (E_s I) \left[ \int_0^\ell w_p w''_{q''} dZ - w_p(\ell) w''_{q''}(\ell) \right]$$
$$
+ \ell(E_s I) w'_p(\ell) w''''(\ell) + K_{\text{Eq-Rub}} w_p(\ell) w_q(\ell) - P(t) \left[ \int_0^\ell w_p w''_q \, dZ - w_p(\ell) w'_q(\ell) \right] \\
- \Omega^2 M_b (r+L/2)^2 w'_p(\ell) w''_q(\ell) \sum_{j=1}^{N_b} \cos^2 (\theta_j + \Omega t) \tag{B53}
$$

$$
[\mathbf{K}_{(U,V)}]_{p,q} = \sum_{i=1}^{N_c} S_{xy}^{b(i)} w_p(z_i) w_q(z_i) + \mu K_{\text{Eq-Rub}} w_p(\ell) w_q(\ell) + 2\psi H\Omega \int_0^\ell w_p w''''_q \, dZ + T(t) \left[ \int_0^\ell w_p w''''_q \, dZ - w'_p(\ell) w''_q(\ell) \right] \\
- \Omega^2 M_b (r+L/2)^2 w'_p(\ell) w''_q(\ell) \sum_{j=1}^{N_b} \sin(\theta_j + \Omega t) \cos(\theta_j + \Omega t) \tag{B54}
$$

In the same manner, the following 2 sub-matrices of Eqs. (B55) and (B56) are repeated j = number of fan blades (N<sub>b</sub>) times:

$$
[\mathbf{K}_{(U,\eta)}^{(j)}]_{p,q} = [\mathbf{K}_{(\eta,U)}^{(j)}]_{q,p} = [-E_s I \cos \beta_0 w_p''''(\ell) \pm \mu K_{\text{Rub}} W_q(L) \cos \beta_L w_p(\ell)] \sin(\theta_j + \Omega t) \\
+ \text{CD} \left[ \begin{aligned}
& W_q'''(0) w_p'(\ell) \text{rsin} \beta_0 \cos(\theta_j + \Omega t) - W_q''(0) w_p'(\ell) \cos \beta_0 \cos(\theta_j + \Omega t) \\
& + W_q'''(0) [w_p(\ell) + (\ell + D_L/2) w_p'(\ell)] \cos \beta_0 \sin(\theta_j + \Omega t)
\end{aligned} \right] \tag{B55}
$$

$$
[\mathbf{K}_{(U,\zeta)}^{(j)}]_{p,q} = [\mathbf{K}_{(\zeta,U)}^{(j)}]_{q,p} = -E_s I w_p'''(\ell) \cos(\theta_j + \Omega t) - \left[ \frac{E_b Ch}{(1-\nu_b^2)} X'_q(0) + K_{\text{Rub}} X_q(L) \right] \left[ w_p(\ell) + (\ell+D_L/2) w'_p(\ell) \right] \cos(\theta_j + \Omega t) \quad (B56)
$$

$$
[\mathbf{K}_{(V,U0)}]_{p,1} = \sum_{i=1}^{N_c} S_{yx}^{b(i)} [w_p(z_i) + z_i w'_p(z_i)] - \mu K_{\text{Eq-Rub}} w_p(\ell) \quad (B57)
$$

$$
[\mathbf{K}_{(V,V0)}]_{p,2} = \sum_{i=1}^{N_c} \left[ \begin{aligned}
& S_{yy}^{b(i)} [w_p(z_i) + z_i w'_p(z_i)] + B_{xx}^{b(i)} w''_p(z_i) \\
& \quad - (E_s I) w'''_p(z_i)
\end{aligned} \right] - (E_s I) w'''_p(\ell) - P(\ell) w'_p(\ell) + K_{\text{Eq-Rub}} w_p(\ell) \quad (B58)
$$

$$
[\mathbf{K}_{(V,U)}]_{p,q} = \sum_{i=1}^{N_c} S_{yx}^{b(i)} w_p(z_i) w_q(z_i) - \mu K_{\text{Eq-Rub}} w_p(\ell) w_q(\ell) - 2\psi H\Omega \int_0^\ell w_p w''''_q \, dZ - T(t) \left[ \int_0^\ell w_p w''''_q \, dZ - w'_p(\ell) w''_q(\ell) \right] \\
- \Omega^2 M_b (r+L/2)^2 w'_p(\ell) w''_q(\ell) \sum_{j=1}^{N_b} \sin(\theta_j + \Omega t) \cos(\theta_j + \Omega t) \quad (B59)
$$

$$
[\mathbf{K}_{(V,V)}]_{p,q} = \sum_{i=1}^{N_c} [S_{yy}^{b(i)} w_p(z_i) w_q(z_i) + B_{xx}^{b(i)} w'_p(z_i) w'_q(z_i)] + (E_s I) \left[ \int_0^\ell w_p w''''_q \, dZ - w_p(\ell) w''''_q(\ell) \right] \\
+ \ell (E_s I) w'_p(\ell) w''''_q(\ell) + K_{\text{Eq-Rub}} w_p(\ell) w_q(\ell) - P(\ell) \left[ \int_0^\ell w_p w''''_q \, dZ - w_p(\ell) w''''_q(\ell) \right] \\
- \Omega^2 M_b (r+L/2)^2 w'_p(\ell) w''''_q(\ell) \sum_{j=1}^{N_b} \sin^2(\theta_j + \Omega t) \quad (B60)
$$

Similarly, for the reasons described earlier, the following sub-matrices of Eqs. (B61)-(B64) are repeated j = number of fan blades (N<sub>b</sub>) times:

$$
[\mathbf{K}_{(V,\eta)}^{(j)}]_{p,q} = [\mathbf{K}_{(\eta,V)}^{(j)}]_{q,p} = [E_s I \cos \beta_0 w_p''''(\ell) \pm \mu K_{\text{Rub}} W_q(L) \cos \beta_L w_p(\ell)] \cos(\theta_j + \Omega t) \\
+ \text{CD} \left[ \begin{aligned}
& W_q'''(0) w_p'(\ell) \text{rsin} \beta_0 \sin (\ell_j + \Omega t) - W_q''(0) w_p'(\ell) \cos \beta_0 \sin (\ell_j + \Omega t) \\
& - W_q'''(0) [w_p(\ell) + (\ell + D_L/2) w_p'(\ell)] \cos \beta_0 \cos (\ell_j + \Omega t)
\end{aligned} \right] \tag{B61}
$$

$$
[\mathbf{K}_{(V,\zeta)}^{(j)}]_{p,q} = [\mathbf{K}_{(\zeta,V)}^{(j)}]_{q,p} = -E_s I w_p'''(\ell) \sin(\theta_j + \Omega t) - \left[ \frac{E_b Ch}{(1-\nu_b^2)} X'_q(0) + K_{\text{Rub}} X_q(L) \right] \\
[w_p(\ell) + (\ell+D_L/2) w'_p(\ell)] \sin(\theta_j + \Omega t) \tag{B62}
$$

$$
[\mathbf{K}_{(\eta,\eta)}^{(j)}]_{p,q} = \sum_{i=1}^{N_c} \frac{z_i}{(\ell+D_L/2)} \left[ \begin{aligned} & S_{xx}^{b(i)} \sin^2 (\theta_j + \Omega t) + S_{yy}^{b(i)} \cos^2 (\theta_j + \Omega t) \\ & - (S_{xy}^{b(i)} + S_{yx}^{b(i)}) \cos (\theta_j + \Omega t) \sin (\theta_j + \Omega t) \end{aligned} \right] \cos^2 \beta_0
$$
$$
\begin{aligned}
& +CD \int_0^L \left[ \begin{aligned} & W_p W_q'''' - 2(\alpha_k^2) W_p W_q'' + (\alpha_k^4) W_p W_q \\ & + (\beta')^2 (3W_p W_q'' - (\alpha_k^2)W_p W_q) \\ & - (\beta')^4 \left(s^2 W_p W_q'' - \frac{(1-v_b)}{2} (\alpha_k^2) s^2 W_p W_q + 2s W_p W_q' \right) \end{aligned} \right] ds \\
& - \frac{(M_b/L)\Omega^2}{2} \int_0^L \left[ \begin{aligned} & (R^2 - r^2 - 2rs - s^2) W_p W_q'' - 2(r+s)W_p W_q' \\ & + 2\cos^2 \beta W_p W_q + (C^2/4) (\alpha_k^2) \sin^2 \beta W_p W_q \end{aligned} \right] ds \\
& +CD \left[ \begin{aligned} & (2 - v_b)(\alpha_k^2) W_p(L) W_q'(L) \\ & - (\beta')^2 [2W_p(L)W_q'(L) - (1 - v_b)(\alpha_k^2)LW_p(L)W_q(L)] \end{aligned} \right] \\
& -CD[v_b(\alpha_k^2) W_p'(L)W_q(L) - (\beta')^2 (W_p'(L)W_q(L) + LW_p'(L)W_q'(L))] \\
& + (F_{\text{Rub}}^{(j)}(t)) \left[ W_p(L)W_q'(L) - \int_0^L W_p(s)W_q''(s)ds \right]
\end{aligned}
\quad (B63)
$$

$$
\begin{aligned}
& [\mathbf{K}_{(\zeta, \zeta')}^{(j)}]_{p,q} = \sum_{i=1}^{N_c} \frac{z_i}{(\ell + D_L/2)} \left[ \begin{aligned}
& S_{xx}^{b(i)} \cos^2 (\theta_j + \Omega t) + S_{yy}^{b(i)} \sin^2 (\theta_j + \Omega t) \\
& + (S_{xy}^{b(i)} + S_{yx}^{b(i)}) \cos (\theta_j + \Omega t) \sin (\theta_j + \Omega t)
\end{aligned} \right] \\
& + [B_{yy}^{b(i)} \cos^2 (\theta_j + \Omega t) + B_{xx}^{b(i)} \sin^2 (\theta_j + \Omega t)] X_p'(0)X_q'(0) \\
& - \frac{E_b Ch}{(1 - v_b^2)} \left[ \int_0^L X_p X_q'' ds - X_p(L)X_q'(L) \right] + K_{\text{Rub}} X_p(L)X_q(L) - (M_b/L)\Omega^2 \int_0^L X_p X_q ds
\end{aligned}
\quad (B64)
$$

In the above equations, viz. Eqs. (B63) and (B64) the contributions of the bearing direct stiffness parameters are included only on the diagonal terms of the sub-matrices $[\mathbf{K}_{(\eta, \eta')}^{(j)}]_{p,q}$ and $[\mathbf{K}_{(\zeta, \zeta')}^{(j)}]_{p,q}$, respectively. Additionally, for any particular jth blade on the overhung fan the magnitude of rub forces $F_{\text{Rub}}^{(j)}(t)$ on the blade tip are computed as follows:

(a) Radial component of blade tip rub-load coming from the abrable layer:

$$[F_{\zeta}^{(j)}]_s = L = -F_{\text{Rub}}^{(j)}(t) = -K_{\text{Rub}}\epsilon_j(t) = -K_{\text{Rub}}[\zeta_0(t) + \zeta(L,t) - \gamma]_j \quad (B65)$$

(b) Tangential component of blade tip rub-load coming from the abrable layer:

$$[F_{\eta}^{(j)}]_s = \mu F_{\text{Rub}}^{(j)}(t) \cos \beta_L = \mu K_{\text{Rub}} \epsilon_j(t) \cos \beta_L = \mu K_{\text{Rub}} [\zeta_0(t) + \zeta(L,t) - \gamma]_j \cos \beta_L \quad (B66)$$

Thus, the column vector related to blade-tip external forces along its longitudinal direction is:

$$\mathbf{p}_{\zeta}^{(j)} = F_{\text{cf}}\mathbf{X}_p(L/2) - F_{\text{Rub}}^{(j)}(t)\mathbf{X}_p(L) \quad (B67)$$

The effect of tip-clearance $\gamma$ is accounted by adding a fictitious force on the right-hand side of the equation which results in,

$$(\mathbf{p}_{\zeta}^{(j)})_{t+\Delta t} = F_{\text{cf}}\mathbf{X}_p(L/2) - F_{\text{Rub}}^{(j)}(t)\mathbf{X}_p(L) + K_{\text{Rub}}\gamma\mathbf{X}_p(L) \quad \text{with } (1 \le p \le n_2) \quad (B68)$$

Similarly, the column vector related to blade-tip external forces in its lateral direction is:

$$(\mathbf{p}_{\eta}^{(j)})_{t+\Delta t} = \mu F_{\text{Rub}}^{(j)}(t) \cos \beta_L \mathbf{W}_p(L) \quad \text{with } (1 \le p \le n_1) \quad (B69)$$

In the above equations, the terms of the series related to the shaft deformation shape functions $w_i(Z)$ and its derivatives along the span of the rotor-shaft ($0 \le Z \le \ell$) are determined as follows:

$$
\left.
\begin{aligned}
w_i(Z) &= [1 - \cos \phi_i Z] \\
w'_i(Z) &= \phi_i \sin(\phi_i Z) \\
w''_i(Z) &= (\phi_i)^2 \cos(\phi_i Z) \\
w'''_i(Z) &= -(\phi_i)^3 \sin(\phi_i Z) \\
w''''_i(Z) &= -(\phi_i)^4 \cos(\phi_i Z)
\end{aligned}
\right\} \quad (B70)
$$

The typical blade deformation vector $\delta_B$ is usually of interest only in rotating coordinates $(\dot{e}_t, \dot{e}_a, \dot{e}_r)$ called tangential-axial-radial system attached with the fan-rotor. The corresponding deformation terms $(\eta, \zeta)$ in the lateral and longitudinal directions of any of the blade thus calculated using Eqs. (B6) and (B7) can easily be converted in the new coordinate system
of interest using,

$$ \hat{\delta}_{\text{B}}(s,c,t) = \eta(s,c,t)\cos(\beta_0 + s\beta')\hat{\mathbf{e}}_t + \eta(s,c,t)\sin(\beta_0 + s\beta')\hat{\mathbf{e}}_a + \zeta(s,t)\hat{\mathbf{e}}_r \quad (B71) $$

References

[1] M.J. Stallone, V. Gallardo, A.F. Storace, L.J. Bach, G. Black, E.F. Gaffney, Blade loss transient dynamic analysis of turbomachinery, AIAA Journal 21 (8) (1983) 1134-1138.

[2] N. Cosme, D. Chevrolet, J. Bonini, B. Peseux, P. Cartraud, Prediction of engine loads and damages due to blade-off event, Proceedings of the 43rd AIAA/ ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Paper No. AIAA-2002-1666, Denver, CO, 22-25 April 2002, 9 pp.

[3] Y.N. Shmotin, D.V. Gabov, A.A. Ryabov, S.S. Kukanov, V.N. Rechkin, Numerical analysis of aircraft engine fan blade-out, Proceedings of the 42nd AIAA/ ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Paper No. AIAA 2066-4620, Sacramento, CA, 9-12 July 2006, 8 pp.

[4] J.B. Husband, Developing an Efficient FEM Structural Simulation of a Fan Blade Off Test in a Turbofan Engine, PhD Thesis, University of Saskatchewan, Saskatoon, Saskatchewan, Canada, October 2007.

[5] S.K. Sinha, S. Dorbala., Dynamic loads in the fan containment structure of a turbofan engine, ASCE Journal of Aerospace Engineering 22 (3) (2009) 260-269.

[6] R. Jain, Prediction of transient loads and perforation of engine casing during blade-off event of fan rotor assembly, Proceedings of the IMPLAST 2010 Conference, Society for Experimental Mechanics, Inc., Providence, RI, October 12-14 2010.

[7] G. von Groll, D.J. Ewins, On the dynamics of windmilling in aero-engines, Proceedings of the IMechE Conference Transactions, Vol. 6, 2000, pp. 721-730.

[8] F.A. Favara, Turbine Engine Continued Rotation and Rotor Locking, Advisory Circular No. 33.74/92-1A, Initiated by ANE-110, US Department of Transportation, Federal Aviation Administration, July 23, 2003.

[9] H.J. Xuan, R.R. Wu, Aeroengine turbine blade containment tests using high-speed rotor spin testing facility, Aerospace Science and Technology 10 (6) (2006) 501-508.

[10] Q. He, H. Xuan, L. Liu, W. Hong, R. Wu, Perforation of aero-engine fan casing by a single rotating blade, Aerospace Science and Technology, 2012, http://dx.doi.org/10.1016/j.ast.2012.01.010 (available online: 14.01.2012).

[11] C. Padova, M.G. Dunn, J. Barton, K. Turner, A. Turner, D. DiTommaso, Casing treatment and blade-tip configuration effects on controlled gas turbine blade tip/shroud rubs at engine conditions, ASME Journal of Turbomachinery 133 (2011) 011016-1-011016-12.

[12] Australian Transport Safety Bureau (ASTB), Examination of Failed Fan Blade Rolls-Royce RB211 Trent 892 Turbofan Engine, Boeing 777-300, A6-EMM, Technical Analysis Report No. 8/01, Occurrence No. 200100445, Reference: BE/200100004 (2001) 21 pp.

[13] O.N. Kirilov, Campbell diagrams of weakly anisotropic flexible rotors, Proceeding of Royal Society A: Mathematical Physical and Engineering Sciences 465 (2009) 2703-2723.

[14] J.S. Rao, R. Sreenivas, Dynamics of asymmetric rotors using solid models, Proceedings of the International Gas Turbine Congress 2003, Paper No. IGTC 2003 Tokyo TS-016, Tokyo, November 2-7 2003, 6 pp.

[15] S.K. Sinha, Stability of viscoelastic rotor-disk system under dynamic axial loads, AIAA Journal 27 (11) (1989) 1653-1655.

[16] S.K. Sinha, On general conditions of rotordynamic stability under combined axial forces and torque, ASME Journal of Applied Mechanics 59 (1992) 225-228.

[17] J. Yoon, J. Lee, Explanation of load torque on rotor stability by internal and external damping, Journal of Aerospace Research from Korea 11 (1993) 481-486.

[18] F.A. Sturla, A. Argento, Free and forced vibrations of a spinning viscoelastic beam, ASME Journal of Vibration and Acoustics 118 (1996) 463-468.

[19] S.H.K. Huang, Coupled shaft torsion and blade bending vibrations of a rotating shaft-blade unit, ASME Journal of Engineering for Gas Turbine and Power 118 (1996) 100-106.

[20] H.C. Sheu, L.W. Chen, Lumped mass model for parametric instability analysis of cantilever shaft-disk systems, Journal of Sound and Vibration 234 (2) (2000) 331-348.

[21] Y.M. Zhang, B.C. Wen, Q.L. Liu, Reliability sensitivity for rotor-stator systems with rubbing, Journal of Sound and Vibration 259 (5) (2003) 1095-1107.

[22] Y.N. Al-Nassar, M. Kalyon, M. Pakdemili, A.O. Al-Bedoor, Stability analysis of rotating blade vibration due to torsional excitation, Journal of Vibration and Control 13 (9-10) (2007) 1379-1391.

[23] Z. Yuan, S. Wang, X. Yue, Z. Li, F. Chu, Dynamic analysis of rotor's radial rub-impact in full degrees of freedom accounting for turborotors's non-linear clearance excitation force, Proceedings of the IMechE, Part C, Journal of Mechanical Engineering Science 222 (9) (2008) 1647-1653.

[24] M. Heidari, D.L. Carlson, S. Sinha, R. Sadeghi, C. Heydari, H. Bayoumi, J. Son, An efficient multi-disciplinary simulation of engine fan-blade out event using MD nastran, Proceedings of the 49th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Paper No. AIAA 2008-2333, Schaumburg, IL, 7-10 April 2008.

[25] C.-T. Chao, T.-L.F. Chiang, C. Chyn, C.-J. Wang, Analysis of vibration in shaft-disc-blades systems due to power faults using a finite-element method, Journal of Chinese Institute of Engineers 34 (2) (2011) 253-264.

[26] M. Legrand, C. Pierre, P. Cartraud, J.-P. Lombard, Two dimensional modeling of a aircraft engine structural bladed disk-casing model interaction, Journal of Sound and Vibration 319 (1-2) (2009) 366-391.

[27] M. Legrand, A. Batailly, C. Pierre, Numerical investigation of abradable coating removal in aircraft engines through plastic constitutive law, ASME Journal of Computational and Nonlinear Dynamics 7 (2012) 011010-1-011010-11.

[28] A. Batailly, M. Legrand, A. Millecamps, F. Garcin, Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact, ASME Journal of Engineering for Gas Turbines and Power 134 (2012) 082504-1-082504-11.

[29] S.-H. Hsieh, J.A. Abel, Comparison of two finite-element approaches for analysis of rotating, bladed-disk assemblies, Journal of Sound and Vibration 182 (1) (1995) 91-107.

[30] G. Genta, On the stability of rotating blade arrays, Journal of Sound and Vibration 273 (2004) 805-836.

[31] S.K. Sinha, Dynamic characteristics of a flexible bladed-rotor with Coulomb damping due to tip-rub, Journal of Sound and Vibration 273 (4-5) (2004) 875-919.

[32] N. Lesaffre, J.-J. Sinou, J.F. Thouverez, Contact analysis of a flexible bladed rotor, European Journal of Mechanics, A/Solids 26 (3) (2007) 541-557.

[33] D. Cao, L. Wang, Y. Chen, W. Huang, Bifurcation and chaos of the bladed overhang rotor system with squeeze film dampers, Science in China, Series E, Technological Sciences 52 (3) (2009) 709-720.

[34] M. Legrand, C. Pierre, B. Peseux, Structural modal interaction of a four degree-of-freedom bladed disk and casing model, ASME Journal of Computational and Nonlinear Dynamics 5 (4) (2010) 041013. (9 pp).

[35] S. Roques, M. Legrand, P. Cartraud, C. Stoisser, C. Pierre, Modeling of a rotor speed transient response with radial rubbing, Journal of Sound and Vibration 329 (5) (2010) 527-546.

[36] M. Gruin, F. Thouverez, L. Blanc, P. Jean, Nonlinear dynamics of a bladed dual-shaft, European Journal of Computational Mechanics 20 (1-4) (2011) 207-225.

[37] Z. Shang, J. Jiang, L. Hong, The global responses characteristics of a rotor/stator rubbing system with dry friction effects, Journal of Sound and Vibration 330 (2011) 2150-2160.

[38] D. Arzina, Vibration Analysis of Compressor Blade Tip-rubbing, Master's Thesis, School of Engineering, Cranfield University, Cranfield, UK, 2011.
[39] R.J. Williams, Simulation of blade casing interaction phenomena in gas turbines resulting from heavy tip rubs using an implicit time marching method, *Proceedings of the ASME Turbo Expo 2011*, ASME Paper No. GT2011-45495, Vancouver, Canada, June 6-10, 2011.
[40] K.E. Turner, M.G. Dunn, C. Padova, Airfoil deflection characteristics during rub events, *ASME Journal of Turbomachinery* 134 (1) (2012) 011018-1-011018-8.
[41] M. Torkhani, L. May, P. Voinis, Light, medium and heavy partial rubs during speed transients of rotating machines: numerical simulation and experimental observation, *Mechanical Systems and Signal Processing* 29 (2012) 45-66.
[42] S. Lahriri, H.I. Weber, I.F. Santos, H. Hartmann, Rotor-stator contact dynamics using a non-ideal drive—theoretical and experimental aspects, *Journal of Sound and Vibration* 331 (20) (2012) 4518-4536.
[43] K.-T. Kim, C.-W. Lee, Dynamic analysis of asymmetric bladed rotors supported by anisotropic stator, *Journal of Sound and Vibration* 331 (24) (2012) 5224-5246.
[44] S.K. Sinha, K.E. Turner, Natural frequencies of pre-twisted blade in a centrifugal force field, *Journal of Sound and Vibration* 330 (11) (2011) 2655-2681.
[45] S.K. Sinha, Non-linear dynamic response of a rotating radial timoshenko beam with periodic pulse loading at the free-end, *International Journal of Non-linear Mechanics* 40 (1) (2005) 113-149.
[46] S.K. Sinha, Combined torsional-bending-axial dynamics of a twisted rotating cantilever beam with contact-impact loads at the free end, *ASME Journal of Applied Mechanics* 74 (2007) 505-522.
[47] K.G. VanDuyn, Bearing Support, United States Patent No. US 7097413 B2, United Technologies Corporation, Hartford, CT, August 29, 2006.
[48] M. Gennaretti, D. Muro, Multiblade reduced-order aerodynamics for state-space aeroelastic modeling of rotors, *AIAA Journal of Aircraft* 49 (2) (2012) 495-502.