# AN ALGORITHM FOR RESPONSE AND STABILITY OF LARGE ORDER NON-LINEAR SYSTEMS - APPLICATION TO ROTOR SYSTEMS

P. Sundararajan, S.T. Noah

► To cite this version:

P. Sundararajan, S.T. Noah. AN ALGORITHM FOR RESPONSE AND STABILITY OF LARGE ORDER NON-LINEAR SYSTEMS - APPLICATION TO ROTOR SYSTEMS. Journal of Sound and Vibration, 1998, 214 (4), 10.1006/jsvi.1998.1614. hal-01353830

HAL Id: hal-01353830

https://hal.science/hal-01353830

Submitted on 14 Aug 2016

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
# AN ALGORITHM FOR RESPONSE AND STABILITY OF LARGE ORDER NON-LINEAR SYSTEMS— APPLICATION TO ROTOR SYSTEMS

P. SUNDARARAJAN AND S. T. NOAH

Mechanical Engineering Department, Texas A&M University, College Station, TX 77843, U.S.A.

A numerical algorithm to calculate the periodic response, stability and bifurcations of a periodically excited non-conservative, Multi-Degree of Freedom (MDOF) system with strong local non-linearities is presented. First, the given large order system is reduced using a fixed-interface component mode synthesis procedure (CMS) in which the degrees of freedom associated with non-linear elements are retained in the physical co-ordinates while all others, whose number far exceeds the number of non-linear DOF, are transformed to modal co-ordinates and reduced using real mode CMS. A shooting and continuation method is then applied to the reduced system to solve for the periodic response. Floquet stability theory is used to calculate stability and bifurcations of the periodic response. The algorithm is applied to study the response to imbalance, stability, and bifurcations of a 24 DOF flexible rotor supported on journal bearings. The results indicate that the proposed algorithm, though approximate, can yield very accurate information about dynamic behavior of large order non-linear systems, even with few numbers of retained component modes. The algorithm, which imposes less demand on computer time and memory, is believed to be of considerable potential in analyzing a variety of practical problems.

## 1. BACKGROUND

The dynamic analysis of multiple shaft rotor-bearing systems can require the assembly and solution of large order sets of ordinary differential equations of motion. Additionally, when one or more non-linear elements such as fluid-film bearing, squeeze-film damper, ball bearing-clearance, fluid seal, etc., are present in the system, which is often the case, the analysis becomes highly complicated. Very often the non-linear equations of motion are linearized in the neighborhood of an operating point and the order of the resulting linear system of equations is reduced using modal transformations and then the reduced model is solved for system response. Though this analysis yields satisfactory results when the rotor system is operating in a weakly non-linear regime (when a fluid-film bearing is operating with low eccentricity (lightly loaded), it is not useful when the rotor-bearing system displays strongly non-linear behavior. For example, fluid-film bearings produce forces which are highly non-linear functions of displacements when operating at high eccentricities (heavily loaded). In such a case, the linearization procedure masks potentially dangerous bifurcations such as subcritical Hopf bifurcations that are predicted by the original non-linear equations of motion. Such a bifurcation which often leads to serious rotordynamics stability problems can only be predicted by retaining higher-order terms. Only a few analysts have dealt with non-linear rotordynamics for large scale systems and even a smaller number of them have dealt with steady state response and stability of such
systems under forced excitations. Very often stability information is predicted using eigenanalysis of the unforced system even when the system is under forced excitation due to imbalance or any other source of excitation. Such a procedure may be valid only for weakly non-linear systems and more general procedures based on modern dynamical systems theory are required to handle strongly non-linear systems. It is the objective of this work to develop numerical procedures to calculate the steady state periodic response, its stability, and bifurcation for a large order nonlinear system under forced excitation.

The large order system of equations which are generated by linear and/or non-linear simulations are usually time consuming to set up and costly to solve in terms of computer time and storage. This is particularly important in design studies which require a completely new system assembly and solution for each new set of system parameters. For most linear transient analyses and for all non-linear analyses, numerical integration is required causing the cost and accuracy of analysis to be strongly dependent on the order of the system equations. Hence, any algorithm for treating large order systems should incorporate some type of reduction scheme to reduce the order of the system.

Few researchers have addressed the dynamic analysis of large order nonlinear system with particular emphasis on the steady state analysis. The method of harmonic balance has been successfully used by Yamauchi [1], Choi and Noah [2], Kim and Noah [3, 4] and Shiau and Jean [5]. Choi and Noah [2] utilized discrete Fourier transform procedures in conjunction with harmonic balance and also included subharmonic response components. Kim and Noah [3, 4] added a dynamic condensation technique to the Harmonic Balance Method to render the problem tractable for large order systems. This condensation procedure reduces the iteration problem so that it involves those co-ordinates which are directly associated with the nonlinear components. The work of Shiau and Jean [5] is similar to that of Choi and Noah [6] and also includes a condensation technique that reduces the problem to co-ordinates that involves only the nonlinear co-ordinates. Their condensation strategy is a generalization of the same strategy developed by McLean and Hahn [7] for analyzing the centered circular orbit response of large rotor systems with squeeze-film dampers.

The collocation method along with a component mode synthesis (CMS) procedure was used by Nataraj and Nelson [8] to calculate the periodic solutions of large order nonlinear systems. The CMS procedure was used to reduce the order of the system and render the problem numerically tractable. Jean and Nelson [9] presented a collocation method that could be utilized for large order nonlinear systems directly in physical co-ordinates and used a reduction procedure similar to Shiau and Jean [5] and McLean and Hahn [10].

Fey *et al.* [11] used a combination of finite-difference and CMS methods to calculate the steady state response of structural systems. The CMS procedure was used to reduce the order of the system, and the finite-difference method with arc-length continuation scheme is applied to calculate the periodic solution and its stability.

Even with reduction procedures, the harmonic balance and collocation methods cited above generate a system of nonlinear algebraic equations whose order is often quite high when more harmonics are included in the assumed response. For example, for every additional harmonic included, the number of unknowns increases by two times the order of the system. Finite difference methods also can result in a large system of nonlinear algebraic equations depending on the number of time steps used. The solution of such large system of equations is time consuming and is prone to convergence problems. The shooting method with an arc-length continuation scheme was successfully used by Sundararajan and Noah [12–14] and Noah and Sundararajan [15] for a low order rotor-bearing system consisting of squeeze-film dampers and journal bearings. The advantage of the shooting method over other methods such as harmonic balance, collocation and finite difference
methods is that the number of resulting nonlinear algebraic equations is never more than the order of the system. The number of resulting nonlinear algebraic equations in the shooting method is not dependent on the number of harmonics present in the response and will be equal to the order of the system. This means if a procedure such as CMS is used to reduce the order of the original system, the number of nonlinear equations from a shooting method will equal the order of the reduced system and never more. Also, from the authors' experience, the nonlinear algebraic equations from the shooting method have better convergence since they are solved by enforcing a two-point boundary value satisfying the periodicity. In periodically forced systems, the shooting procedure can be started with an assumed period equal to the period of excitation. Since shooting methods require calculation of a Jacobian at every time increment, the Jacobian matrix upon convergence can be used to form the Floquet matrix which in turn will yield local stability and bifurcation information. Because of these advantages, the shooting method is proposed to be used for this study to calculate the periodic solution. A pseudo-arc length continuation scheme is embedded into the shooting scheme to calculate unstable periodic solution branches involving saddle-node bifurcations [16].

For large order systems, Crandall and Yeh [17] proposed a method of automatically generating modes based on Raleigh-Ritz procedure and then reducing the order of the system. This procedure did not uncouple the equations of motions. Glasgow and Nelson [18] were one of the first to recognize the advantages of CMS procedures for reduction of large order rotordynamic systems and used it to calculate the eigenvalues and stability of a rotor system. Nelson et al. [19] further extended the CMS procedure to calculate the transient response of a nonlinear rotor system. The main advantages of the CMS procedure (Childs [20]) in addition to substantial reduction in the total number of degrees of freedom are (1) the CMS procedure lends itself to a transient analysis which is very important for the operation and analysis of many machinery, (2) it can be directly applied to multiple shafts with multiple interconnections, (3) it can be adapted easily to general finite-element models for machinery rotor/housing. Because of its generality and versatility, a CMS procedure is proposed to be used in this work to reduce the order of the system. A brief review of component mode synthesis procedures is given below.

Summaries of the various developments in component mode synthesis were presented by Hou [21], Craig and Chung [22], Nelson et al. [19] and Craig [23]. The initial development of the method of component mode synthesis is due to Hurty [24] with subsequent simplification of the formulation by Craig and Bampton [25]. Glasgow and Nelson [18] developed a CMS approach for rotordynamic systems which utilizes component constraint modes as well as constrained precessional modes, thereby allowing a full modal transformation for general second order systems analogous to the approach by Craig and Bampton [25] for lightly damped symmetric systems. The fixed-interface component mode method used in their work requires that all of the original interface displacement co-ordinates be retained in the final coupled structure model but produces a very accurate eigen solutions for a given number of degrees of freedom [23]. Free-interface component mode synthesis methods are attractive when a structure is assembled using component modes obtained from modal testing. Though these methods have been successfully used by a few authors [26, 27], residual flexibility and inertia-relief modes are to be included to achieve comparable accuracy and convergence as that of the fixed-interface component methods [28, 29] and this will necessitate considerable additional computer programming and book keeping for large order systems. Since fixed-interface CMS procedures have been known to be accurate, it is used in the current work.

The procedure developed by Glasgow and Nelson [18] utilizes fixed component modes for the rotating assembly that includes gyroscopic effects. A primary disadvantage of this
method is that the modes are speed dependent, due to the gyroscopics, and require modification at each spin speed. Also, the equations of motion involve physical co-ordinates which are real values and modal co-ordinates which are complex, and this requires a numerical integration of a system of complex variables. The use of complex modes is justified by stating that the modes are true component modes of the rotating system so that only a few modes need to be retained at any speed in order to obtain a highly accurate model. One of the objectives of this work is to compare the effectiveness of real (planar) modes with that of complex modes in predicting the modal characteristics of a rotor bearing system.

A study using combination of the shooting/continuation method with a CMS procedure has not been reported previously. Because of the cited advantages of both methods, a combination of these methods for engineering application is attractive and is hence the subject of this study. A detailed mathematical formulation of shooting and continuation methods for non-autonomous systems are given by Keller [16], Seydel [30], and Parker and Chua [31]. A detailed exposition of CMS procedures for rotor-bearing systems can be found in Glasgow and Nelson [18], Nelson *et al.* [19] and Childs [20]. Hence, the details of these methods will not be repeated here. Only relevant details with regard to the current application will be given here. In the following section a general mathematical development of the combined procedure involving CMS, shooting, and continuation methods is presented.

## 2. MATHEMATICAL DEVELOPMENT

The proposed combined algorithm first uses a fixed-interface CMS procedure to substantially reduce the order of the system and then applies shooting and continuation schemes to calculate the periodic solution, its stability, and bifurcations. First a finite-element model is developed for the rotor system and the degrees of freedom associated with nonlinear supports are identified. These degrees of freedom will be called boundary co-ordinates and will be retained in the analysis in the original physical co-ordinates themselves. The remaining degrees of freedom not associated with the nonlinear elements will be called internal co-ordinates. These degrees of freedom will be converted to modal co-ordinates and co-ordinates corresponding to higher modes will be subsequently discarded thereby effectively reducing the total number of degrees of freedom. Since in a rotor system, the forcing function is predominantly at the running speed (due to imbalance), it is often enough to retain only modes with frequencies up to and including the running speed. Stability and bifurcations of the periodic response for the reduced system is then calculated using the shooting and continuation scheme. The shooting and continuation schemes generate a matrix called Monodromy matrix (also called Jacobean) and the eigenvalues of this matrix help in determining the stability of the periodic response. If all the eigenvalues are within the unit circle, then the periodic solution is stable and if any one of them is outside the circle, then the solution is unstable due to a bifurcation. The way through which an eigenvalue or a pair of eigenvalues cross the unit circle determines the nature of instability. A reader might refer to a number of excellent texts available on this subject such as Balachandran and Nayfeh [32].

The following steps describe the mathematical details of system order reduction using the CMS procedure.

1. Symmetric mass, damping, and stiffness matrices for the system are assembled using a finite-element formulation or any other method. For reasons to be explained later, the non-symmetric stiffness and damping coefficients from elements such as bearings, and seals and gyroscopic forces are not included in the assembly procedure at this time.
(2) The boundary co-ordinates ($x_b$) are identified. The boundary co-ordinates are the degrees of freedom associated with nonlinear elements. The remaining degrees of freedom are called internal co-ordinates ($x_i$).

(3) The mass, stiffness, damping matrices and the force vectors are rearranged using simple transformations to determine submatrices of each of the above matrices/vectors associated with boundary and interior co-ordinates, respectively:

$$ \begin{bmatrix} M_{bb} & M_{bi} \\ M_{ib} & M_{ii} \end{bmatrix} \begin{Bmatrix} \ddot{x}_b \\ \ddot{x}_i \end{Bmatrix} + \begin{bmatrix} C_{bb} & C_{bi} \\ C_{ib} & C_{ii} \end{bmatrix} \begin{Bmatrix} \dot{x}_b \\ \dot{x}_i \end{Bmatrix} + \begin{bmatrix} K_{bb} & K_{bi} \\ K_{ib} & K_{ii} \end{bmatrix} \begin{Bmatrix} x_b \\ x_i \end{Bmatrix} = \begin{Bmatrix} F_b \\ F_i \end{Bmatrix} + \begin{Bmatrix} F_{nsb} \\ F_{nsi} \end{Bmatrix}. \quad (1) $$

The subscript $s$ indicates that matrices are symmetric. $F_b$ denotes forces acting at the boundary co-ordinates which can be nonlinear functions of boundary co-ordinates. $F_i$ denotes forces acting at the interior co-ordinates and can be unbalance forces or any other excitation. $F_{nsb}$ and $F_{nsi}$ are the forces from elements such as a bearing or a seal with non-symmetric stiffness and damping properties or from gyroscopics. The gyroscopic matrix is defined as

$$ \begin{bmatrix} 0 & \omega[G] \\ -\omega[G] & 0 \end{bmatrix} \begin{Bmatrix} \dot{u}_x \\ \dot{u}_y \end{Bmatrix} $$

where $[G]$ defines the gyroscopic coupling properties, $\dot{u}_x$, $\dot{u}_y$ are velocities of degrees of freedom associated with disk nodes.

(4) The main idea of the fixed-interface component mode synthesis (CMS) procedure is to express the displacement of any point in a component as the superposition of two types of displacement modes, that are (i) constrained normal modes—displacement relative to the fixed component boundaries, and (ii) static modes—displacement produced by displacing only the boundary co-ordinates. The constrained normal modes are obtained by setting the boundary co-ordinate state vector in equation (1) to zero and then obtaining the free-vibration problem for the interior co-ordinates:

$$ M_{ii}\ddot{x}_i + C_{ii}\dot{x}_i + K_{ii}x_i = 0. \quad (2) $$

The eigenvalues and eigenvectors of this system are complex in general. The complex eigenvectors, as seen from the works of Glasgow and Nelson [18] and Nelson et al. [19] when used for modal truncation, results in modal co-ordinates that are complex and require extra computer memory and time in addition to requiring numerical integration of complex variables (transient analysis). An alternative is not to include matrix $C_{ii}$, at the expense of a small loss of accuracy, while calculating the eigenvectors of the above system in which case the resulting eigenvectors are real. This is because matrix $M_{ii}$ is positive-definite and $K_{ii}$ is symmetric in general. The main advantage of using real modes, despite being approximate, is that they are not a function of speed and hence need to be calculated only once. On the contrary, if the non-symmetric stiffness and damping properties of a seal or a bearing in the original system were included in the calculation of the eigenvectors required for the CMS procedure, then the eigenvectors would have to be calculated at every speed since mass, stiffness and damping properties of fluid-film bearings, seals etc. are functions of speed. The assembly procedure has to be repeated at every speed. Also, the procedure yields complex eigenvalues and eigenvectors. Because of the simplicity of real modes, it is proposed to examine their usefulness in calculating unbalance response, stability, and bifurcations of a large order nonlinear system under periodic excitation.
Let the real eigenvectors of the system,

$$M_{ii}\ddot{x}_i + K_{ii}x_i = 0, \quad (3)$$

form an eigenvector matrix denoted by **A**.

(5) The static vectors which are displacements produced by displacing the boundary co-ordinates are given by the matrix **B**,

$$x_i = Bx_b, \quad B = -K_{ii}^{-1}K_{ib}. \quad (4)$$

(6) Now the desired co-ordinate transformation mentioned in step 4 is achieved from, $x_i = Bx_b + Aq$, where **q** is the modal co-ordinate vector.

(7) The proposed transformation does not change the boundary co-ordinates; hence, the complete transformation is

$$\begin{Bmatrix} \mathbf{x}_b \\ \mathbf{x}_i \end{Bmatrix} = \begin{bmatrix} \mathbf{I}_{bb} & \mathbf{0}_{bi} \\ \mathbf{B}_{ib} & \mathbf{A}_{ii} \end{bmatrix} \begin{Bmatrix} \mathbf{x}_b \\ \mathbf{q}_i \end{Bmatrix}, \quad \begin{Bmatrix} \mathbf{x}_b \\ \mathbf{x}_i \end{Bmatrix} = \mathbf{D} \begin{Bmatrix} \mathbf{x}_b \\ \mathbf{q}_i \end{Bmatrix}, \quad (5)$$

where $\mathbf{I}_{bb}$ is a unit matrix. **D** takes its meaning from equation (5). Substituting in equation (1) and multiplying the resulting equation by $\mathbf{D}^T$ yields

$$\begin{bmatrix} \mathbf{M}_{11} & \mathbf{M}_{12} \\ \mathbf{M}_{21} & \mathbf{M}_{22} \end{bmatrix} \begin{Bmatrix} \ddot{\mathbf{x}}_b \\ \ddot{\mathbf{q}} \end{Bmatrix} + \begin{bmatrix} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{bmatrix} \begin{Bmatrix} \dot{\mathbf{x}}_b \\ \dot{\mathbf{q}} \end{Bmatrix} + \begin{bmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} \\ \mathbf{K}_{21} & \mathbf{K}_{22} \end{bmatrix} \begin{Bmatrix} \mathbf{x}_b \\ \mathbf{q} \end{Bmatrix} = \begin{Bmatrix} \mathbf{F1} \\ \mathbf{F2} \end{Bmatrix}. \quad (6)$$

By choosing only the modes that are of interest in matrix **A**, one can drastically reduce the size of **D** and hence that of the system. Typically, only modes whose natural frequencies are within or moderately above the running speed of the rotor are retained. While the above equations are coupled, they are much smaller in number compared to the original system in equation (1). The reduced model can be used for stability, synchronous-response, or transient analysis.

(8) The reduced system is to be solved using the shooting and continuation schemes. First, the system has to be written in first order form for ease of numerical integration. This yields

$$\begin{gathered} \begin{bmatrix} \mathbf{0} & \mathbf{M} \\ \mathbf{M} & \mathbf{C} \end{bmatrix} \dot{\mathbf{h}} + \begin{bmatrix} -\mathbf{M} & \mathbf{0} \\ \mathbf{0} & \mathbf{K} \end{bmatrix} \mathbf{h} = \begin{Bmatrix} \mathbf{0} \\ \mathbf{F} \end{Bmatrix}, \\ \dot{\mathbf{x}}_b = \begin{bmatrix} \dot{\mathbf{q}} \\ \mathbf{x}_b \\ \mathbf{q} \end{bmatrix} \\ A_1 \dot{\mathbf{h}} + B_1 \mathbf{h} = F_1, \end{gathered} \quad (7)$$

where matrices $\mathbf{A}_1$, $\mathbf{B}_1$, and $\mathbf{F}_1$ take their meaning from equation (7).

(9) It is convenient in this first order development to introduce the transformation

$$\mathbf{h} = \mathbf{a}\mathbf{y}, \quad \mathbf{y} = (\dot{\mathbf{x}}_b, \mathbf{x}_b, \dot{\mathbf{q}}, \mathbf{q})^T. \quad (8)$$

The above transformation groups the boundary degrees of freedom together and hence also the right-hand side for the nonlinear forces at the boundary co-ordinates which makes the Jacobian formulation easier.
The above substitution yields a system of first order equations,

$$E\dot{y} + Fy = Q, \qquad (9)$$

where **E**, **F**, **Q** are corresponding transformed matrices.

The system is now in a form amenable to the shooting and continuation schemes briefly described below. The shooting and continuation methods are used to solve for the steady state periodic response, its stability, and bifurcations when a parameter (speed) is varied.

### 3. SHOOTING METHOD FOR PERIODIC SOLUTION OF PERIODICALLY FORCED SYSTEMS

A shooting method to locate a periodic solution of a nonlinear system subjected to external excitation (non-autonomous system) is presented in this section. The following exposition is a modification of the procedure of Aluko and Chang [33] for autonomous systems and closely follows that in the paper by Sundararajan and Noah [34]. Periodic solutions are sought for a system described by a given system of differential equations stated in the first order form:

$$\begin{gathered} \dot{\mathbf{x}}(\tau) = \mathbf{f}(\tau, \mathbf{x}(\tau), \lambda), \\ \mathbf{x}, \mathbf{f} \in \mathbb{R}^n; \mathbf{x} = (x_1, x_2, \dots, x_n), \\ \lambda, \tau \in \mathbb{R}, \tau \ge 0, \end{gathered} \qquad (10)$$

where the ‘over dot’ signifies differentiation with respect to time $\tau$. The function $\mathbf{f}$ is an explicit function of time indicating that the system considered is non-autonomous. The dependence of the periodic solutions of equation (10) on the bifurcation parameter $\lambda$ will be investigated.

Any periodic solution of equation (10) must satisfy

$$\mathbf{x}(0) = \mathbf{x}(T), \qquad (11)$$

where $T$ is the minimum period of the response. It is often convenient to normalize the period to unity so that the integration time interval is $[0, 1]$. With this, the differential equation (10) and its imposed condition, equation (11), are transformed to

$$\dot{\mathbf{x}}(t) = T\mathbf{f}(t, \mathbf{x}(t), \lambda) \equiv \mathbf{g}(t, \mathbf{x}: \lambda), \quad \mathbf{x}(0) = \mathbf{x}(1), \qquad (12)$$

where $\mathbf{g} = [g_1, g_2, \dots, g_n]^T$, $t = \tau/T$.

The original problem has now been replaced by a boundary value problem with non-separated boundary conditions, which will prove to be useful in the continuation procedure to be explained later.

#### 3.1. NEWTON-FOX PROCEDURE

Keller's outline [16] for the solution of a generalized two-point boundary problem, adopted to the present case for locating periodic solutions [33], is used here. The boundary value problem can be restated as

$$\dot{\mathbf{x}}(t) = \mathbf{g}(t, \mathbf{x}: \lambda), \quad \mathbf{E}\mathbf{x}(0) + \mathbf{F}\mathbf{x}(1) = 0, \qquad (13a, b)$$

where **E** is a $n \times n$ unit matrix and $\mathbf{F} = -\mathbf{E}$. Consider an initial value problem related to equation (13a),

$$\dot{\mathbf{u}} = \mathbf{g}(t, \mathbf{u}: \lambda), \quad \mathbf{u}(0) = \mathbf{s}, \quad \mathbf{u}, \mathbf{s} \in \mathbb{R}^n. \qquad (14)$$
Denoting the solution of equation (14) by $\mathbf{u}(\mathbf{s}, \lambda, t)$, we seek $\mathbf{s}$ such that

$$ \phi(\mathbf{s}, \lambda) = 0, \qquad (15) $$

where the boundary conditions residual function $\phi$ is defined as

$$ \phi(\mathbf{s}, \lambda) = \mathbf{E}\mathbf{u}(\mathbf{s}, \lambda, 0) + \mathbf{F}\mathbf{u}(\mathbf{s}, \lambda, 1). \qquad (16) $$

Obviously, solution $\mathbf{s}$ will be an $n$-dimensional vector on the periodic solution with the assumed period $T$. In case such a solution cannot be found, the algorithm has to be repeated for increasing integer multiples of $T$ until a solution is found. With a given period $T$ and parameter $\lambda$, there are $n$ unknowns (components of vector $\mathbf{s}$) in as many equations.

The basic Newton-Fox shooting algorithm may be outlined as follows:

$\mathbf{s}^0$—initial guess of the solution, $\mathbf{s}^{n+1} = \mathbf{s}^n + \Delta\mathbf{s}^n$,

where the corrections $\Delta\mathbf{s}^n$ solve the linear system of equations viz.,

$$ \mathbf{Q}(\mathbf{s}^n, \lambda_0)(\Delta\mathbf{s}^n) = -\phi(\mathbf{s}^n, \lambda_0). \qquad (17) $$

The $n$-order square matrix $\mathbf{Q}$ in equation (17) is the Jacobian of $\phi(\mathbf{s}, \lambda)$ with respect to $\mathbf{s}$. This can be written as

$$ \mathbf{Q}(\mathbf{s}, \lambda_0) = \left[ \frac{\partial \phi}{\partial \mathbf{s}} \right] = [\mathbf{E} + \mathbf{F}\mathbf{W}(\mathbf{s}, \lambda, 1)], \qquad (18) $$

where the $n \times n$ matrix $\mathbf{W}(\mathbf{s}, \lambda, t)$ is obtained by solving the following variational set of equations associated with equation (14):

$$ \dot{\mathbf{W}} = \mathbf{A}(t; \mathbf{s}, \lambda)\mathbf{W}, \quad \mathbf{W}(\mathbf{s}, \lambda, 0) = \mathbf{I}. \qquad (19) $$

The $n \times n$ matrix $\mathbf{A}$ is defined as

$$ \mathbf{A}(t; \mathbf{s}, \lambda) = \frac{\partial \mathbf{g}}{\partial \mathbf{u}}(t; \mathbf{u}(\mathbf{s}, \lambda, t), \lambda). \qquad (20) $$

The $n$-order matrix $\mathbf{A}$ is thus the Jacobian of function $\mathbf{g}$ with respect to $\mathbf{u}$ evaluated along the periodic solution.

The above iterative procedure for computation of a point on the periodic solution at $\lambda = \lambda_0$ can be summarized as follows.

For a given $\mathbf{s}^n$, we simultaneously integrate the given equations of motion (14), and the variational set (19) to get $\mathbf{W}(\mathbf{s}, \lambda, 1)$ and $\phi(\mathbf{s}^n, \lambda_0)$ and solve the linear system (17) for $\mathbf{s}^n$ at each iteration. When $\Delta\mathbf{s}^n$ is sufficiently low, convergence to the periodic solution is achieved and the iterations are terminated. In the present work the following relative error condition [16] was used to check convergence (symbol $\langle x, y \rangle$ denotes inner-product of $x$ and $y$),

$$ \frac{\langle \Delta \mathbf{s}^n, \Delta \mathbf{s}^n \rangle}{\langle \mathbf{s}^n, \mathbf{s}^n \rangle} \le 1.0E-9. \qquad (21) $$

## 4. CONTINUATION OF A PERIODIC SOLUTIONS—PSEUDO-ARC LENGTH ALGORITHM

The problem of continuation of solutions in general is that of computing whole solution branches, i.e., obtaining a solution at $\lambda = \lambda_{j+1}$ given that the accurate solution at $\lambda = \lambda_j$ is known. Several continuation schemes have been outlined by Seydel [30]. Among them
the pseudo-arc length continuation scheme has been widely used by several authors [35, 36] due to its accuracy and simplicity. To this end, Keller's algorithm [16] for algebraic continuation has been adopted in the present work for determining the periodic solutions. In order to use this algorithm, the problem of continuation of the periodic solution is converted to the problem of continuation of the fixed points of the Poincare map of the system. In other words, the continuation may be considered to be one of the solutions of the previously formulated system of algebraic equations (15).

The basic idea in pseudo-arc length continuation is to drop the natural parameterization by $\lambda$ and use some other parameterization. Consider the equation

$$ \phi(\mathbf{s}, \lambda) = 0, $$

where $\phi$ is the boundary condition residual function defined in equation (15). If $(\mathbf{s}^0, \lambda^0)$ is any point on a regular path and $(\mathbf{s}^0, \mathbf{0})$ is the unit tangent to the path, the scalar normalization is adjoined to equation (15) so that

$$ N(s, \lambda, \Delta L) = \dot{\mathbf{s}}^0(\mathbf{s} - \mathbf{s}^0) + \dot{\lambda}^0(\lambda - \lambda^0) - \Delta L = 0. \quad (22) $$

This is the equation of plane, which is perpendicular to the tangent $(\mathbf{s}^0, \mathbf{0})$ at a distance $\Delta L$ from $(\mathbf{s}^0, \lambda^0)$. This plane will intersect the curve $\Gamma(L)$, if $\Delta L$ and the curvature of $\Gamma$ are not too large. That is, equations (15) and (22) are solved simultaneously for $(\mathbf{s}(L), \lambda(L))$. Using Newton's method, this leads to a linear system:

$$ \begin{bmatrix} \phi_s^v & \phi_\lambda^v \\ \dot{\mathbf{s}}^0 & \dot{\lambda}^0 \end{bmatrix} \begin{bmatrix} \Delta \mathbf{s}^v \\ \Delta \lambda^v \end{bmatrix} = - \begin{bmatrix} \phi^v \\ N^v \end{bmatrix}. \quad (23) $$

Here, $\phi_s^v = \phi_s^v(\mathbf{s}^v, \lambda^v)$, $\phi_\lambda^v = \phi_\lambda^v(\mathbf{s}^v, \lambda^v)$, and the iterates are $\mathbf{s}^{v+1} = \mathbf{s}^v + \Delta\mathbf{s}^v$ and $\lambda^{v+1} = \lambda^v + \Delta\lambda^v$. The coefficient matrix in the above relation is always non-singular. However, whenever $\phi_s$ is singular or almost singular a bordering algorithm [16] can be used.

It may be noticed that the continuation routine requires the evaluation of $\phi_\lambda$ at each iteration. The formulation is

$$ \frac{\partial \phi}{\partial \lambda} = \mathbf{FZ}(1), \quad (24) $$

where $\mathbf{F}$ is the same boundary condition matrix defined in equation (13b) and $\mathbf{Z}(t) = \partial \mathbf{u}/\partial \lambda$ is the $n$-vector solution of the variational system,

$$ \dot{\mathbf{Z}} = \mathbf{A}\mathbf{Z} - \mathbf{R}, \quad \mathbf{Z}(\mathbf{0}) = \mathbf{0}. \quad (25) $$

The $n \times n$ matrix $\mathbf{A}$ already was defined in equation (20) while the $n$-vector $\mathbf{R}$ is

$$ \mathbf{R} = \frac{\partial \mathbf{g}}{\partial \lambda}. \quad (26) $$

Thus, when expressions for $\mathbf{R}$ are available, they can be inserted into equation (19) and integrated with zero initial conditions. The product of the constant matrix $\mathbf{F}$ and the matrix $\mathbf{Z}(t)$ when computed at the end of the period then gives $\phi_\lambda$ in equation (23).

Determination of the unit tangent $(\mathbf{s})$ is done only up to its sign. The choice of the sign is crucial in order to avoid going backwards or getting trapped between two points. If the tangent at the previous solution is $(\mathbf{s}^0)$, one chooses the sign for $(\mathbf{s})$ so as to make the inner product $\langle(\mathbf{s}^0, \mathbf{0}), (\mathbf{s})\rangle$ positive. This unit tangent vectors $(\mathbf{s}^0, \mathbf{0})$ must satisfy

$$ \phi_s^0 \dot{\mathbf{s}}^0 + \phi_\lambda^0 \dot{\lambda}^0 = 0 \quad (27a) $$
and

$$|\dot{\mathbf{u}}^0|^2 + |\dot{\lambda}^0|^2 = 1. \qquad (27b)$$

Let us consider regular points where $\phi_s^0$ is non-singular. The $\phi_0$ is found from

$$\phi_s^0 \gamma_0 = -\phi_\lambda^0. \qquad (28)$$

Then set,

$$\dot{\mathbf{s}}^0 = a\gamma_0, \qquad \dot{\lambda}_0 = a, \qquad (29)$$

where $\mathbf{a}$ is determined from equation (30) as

$$a = \pm \frac{1}{\sqrt{1 + |\gamma_0|^2}}. \qquad (30)$$

The sign of $\mathbf{a}$ is chosen so that the orientation of the path is preserved, if $(s_{-1}, \lambda_{-1})$ is the preceding tangent vector, then we require

$$a[\dot{\mathbf{s}}_{-1}^T \gamma_0 + \dot{\lambda}_{-1}] > 0. \qquad (31)$$

In the sections that follow, first numerical examples to compare the accuracy of eigenvalues predicted by the CMS procedure using complex modes with those produced by the CMS procedure using real modes are presented for MDOF rotor-bearing systems. The objective here is to see if real modes can produce reasonably accurate results compared to those obtained using complex modes. If the real modes are accurate enough, then they can be used for obtaining steady state response and stability of the system as well. Use of real modes over complex modes could save computer time and memory, and programming complexity. Numerical examples using real modes for obtaining the steady state response and stability are presented in a subsequent section.

## 5. NUMERICAL EXAMPLES—EIGENVALUE ANALYSIS

Two example rotor-bearing systems are used to compare the accuracy of eigenvalues predicted by the CMS procedure using real modes with those calculated using complex modes. The first example is a single uniform shaft used by Glasgow and Nelson [18]. The shaft is a 10·27 cm(4") diameter steel member, 127 cm (50") long, supported on identical bearings. This example was chosen so that the results using the proposed algorithm can be bench marked with those published in the above references. The shaft is modelled as an assemblage of five equal length finite elements. Whirl modes and stability results are presented for this simple model. The rotor model has 24 displacement degrees of freedom, 4 of which are designated as boundary co-ordinates. Thus, there are 20 constrained precessional modes available for truncation. For example, with 16 constrained precessional modes truncated, the original 24 degrees of freedom system is then approximated by an 8 degree of freedom model. In Case 1, a simple shaft supported on isotropic bearings is analyzed. In Case 2 the same beam is used but in an overhung configuration with additional masses and rotational inertias (disks), and stiffness, damping, and mass coefficients from a liquid annular seal. Case 2 thus as pronounced gyroscopic and external damping forces in the model. Numerical values of the physical and support parameters are given below for both cases. Truncated whirl frequency analyses were performed for these cases.
Case 1 (refer to Figure 1): undamped isotropic bearings with stiffness $K = 1.753 \times 10^7$ N/m (1.0 $10^5$ lb/in.) are located at the two ends of a uniform shaft [18]. The shaft is a 10.27 cm (4") diameter steel member with a length of 127 cm (50"). The shaft is modelled using five finite elements with six nodes.

Case 2 (refer to Figure 2): same rotor as in Case (1) is used with the following changes in the configuration. (1) The isotropic bearings are located at nodes 1 and 2 where node 1 is at one end of the rotor and node 2 is 25 cm (10") away from it. (2) A liquid annular seal with following characteristics is located on node 3. These characteristics were taken from a numerical example presented by Childs [20]. $K = 1.567 \times 10^7$ N/m (89 542 lb/in.), $k = 0.473 \times 10^7$ N/m (27 028 lb/in.), $C = 0.3 \times 10^5$ Ns/m (171.42 lb-s/in.), $c = 0.210 \times 10^4$ N/m (12 lb/in.), $M = 6.68$ kg.

Two similar disks with following mass properties are added at nodes 5 and 6: $M = 4.83$ kg (0.028 slug), $I_p = 2.4 \times 10^{-3}$ kg m$^2$ (0.224 slug in$^2$), $I_T = 1.23 \times 10^{-3}$ kg m$^2$ (0.114 slug in$^2$). The whirl frequency occurring at a rotation speed of 4000 r.p.m. for various levels of mode truncation are tabulated in Tables 1 and 2 for both Cases 1 and 2.

Table 1 shows the calculated eigenvalues using real and complex modes are shown opposite to the values reported by Glasgow and Nelson [8] who used complex modes for the uniform beam rotor supported on isotropic bearings (Case 1). It can be seen that the eigenvalues predicted using real and complex modes are nearly identical for this case, irrespective of the number of modes retained. The difference is less than 1% for all cases. The calculated eigenvalues compare very well with Glasgow and Nelson [18] and are in fact slightly more accurate. This is probably due to better numerical accuracies attainable in current day computers. The Case 1 study provided a benchmark since the eigenvalues were verified to be in agreement with those in Glasgow and Nelson [18].

Predicted eigenvalues for Case 2 using real and complex modes are given in Table 2. The rotor in this case has significant gyroscopic effects and external direct/cross-coupled damping, stiffness, and mass coefficients that are due to a liquid seal. The results show that when the number of precessional modes retained in more than 4, the real and complex mode CMS procedures yield almost identical eigenvalues. When the number of retained modes is 4 or less, the eigenvalues corresponding to higher modes become more and more inaccurate. It is also seen that complex modes produce slightly more accurate values of the real part of the eigenvalues (modal damping) for higher modes. However, the first four eigenvalues are still quite accurate within 5% of the values obtained using a higher number of retained modes by both real and complex CMS procedures.

Few more cases of rotor supported on bearings with unsymmetric dynamic coefficients were also studied. The results tentatively show that when gyroscopic and unsymmetric properties of rotor bearing system are negligible, the real and complex eigenvectors produce almost identical results and both are equally accurate. When the number of retained modes is less, both the procedures produce equally inaccurate results for higher modes. For Case 1, a maximum inaccuracy of 5% was observed for the lowest number of retained modes. When gyroscopic and external damping effects are significant, both real and complex modes produce accurate results when the number of retained precessional modes is high. When the number of retained modes is less, only lower modes are predicted accurately. Both produce almost equally inaccurate results for higher modes although complex modes produce more accurate values of the real part of the eigenvalue (modal damping). The author believes complex modes will yield more accurate results when the excitation frequencies are higher than the first few eigenvalues. Since such high excitation frequencies are rare in a rotating machinery; real modes CMS should be sufficient for calculation of eigenvalues, transient analysis and steady state analysis of most systems. It is recommended that real modes be used whenever possible and complex modes are only
<table><caption>TABLE 1
Comparison of eigenvalues from real and complex CMS procedures—Case 1</caption><thead><tr><th rowspan="3">Mode</th><th colspan="12">Number of constrained precessional modes retained</th></tr><tr><th colspan="3">20</th><th colspan="3">12</th><th colspan="3">8</th><th colspan="3">4</th><th colspan="3">2</th></tr><tr><th>Real</th><th>Com</th><th>Nel</th><th>Real</th><th>Com</th><th>Nel</th><th>Real</th><th>Com</th><th>Nel</th><th>Real</th><th>Com</th><th>Nel</th><th>Real</th><th>Com</th><th>Nel</th></tr></thead><tbody><tr><td>1B</td><td>520.6</td><td>520.6</td><td>520.5</td><td>520.6</td><td>520.6</td><td>518.2</td><td>520.7</td><td>520.6</td><td>520.7</td><td>520.7</td><td>520.4</td><td>520.7</td><td>520.7</td><td>522.4</td></tr><tr><td>1F</td><td>521.3</td><td>521.3</td><td>521.0</td><td>521.3</td><td>521.3</td><td>520.0</td><td>521.3</td><td>521.3</td><td>521.2</td><td>521.3</td><td>521.4</td><td>523.0</td><td>521.4</td><td>523.0</td></tr><tr><td>2B</td><td>1093</td><td>1093</td><td>1093</td><td>1093</td><td>1090</td><td>1093</td><td>1093</td><td>1091</td><td>1094</td><td>1094</td><td>1095</td><td>1137</td><td>1137</td><td>1137</td></tr><tr><td>2F</td><td>1097</td><td>1097</td><td>1097</td><td>1097</td><td>1094</td><td>1097</td><td>1098</td><td>1094</td><td>1097</td><td>1098</td><td>1098</td><td>1141</td><td>1141</td><td>1141</td></tr><tr><td>3B</td><td>2237</td><td>2237</td><td>2238</td><td>2238</td><td>2311</td><td>2240</td><td>2240</td><td>2221</td><td>2221</td><td>2290</td><td>2290</td><td>2144</td><td>2290</td><td>2144</td></tr><tr><td>3F</td><td>2255</td><td>2255</td><td>2253</td><td>2255</td><td>2255</td><td>2327</td><td>2257</td><td>2257</td><td>2236</td><td>2309</td><td>2309</td><td>2156</td><td>2309</td><td>2159</td></tr><tr><td>4B</td><td>5058</td><td>5058</td><td>5058</td><td>5059</td><td>5059</td><td>5384</td><td>5067</td><td>5067</td><td>5168</td><td>5259</td><td>5168</td><td>5048</td><td>—</td><td>—</td></tr><tr><td>4F</td><td>5098</td><td>5058</td><td>5093</td><td>5099</td><td>5099</td><td>5424</td><td>5107</td><td>5107</td><td>5297</td><td>5211</td><td>5211</td><td>5085</td><td>—</td><td>—</td></tr></tbody></table>

Notes: Real = real modes, Com = complex modes, Nel = Glasgow and Nelson [18]; Eigenvalues in Table 1 are expressed in rad/s, B, F = backward and forward modes.
<table><caption>TABLE 2
Comparison of eigenvalues from real and complex CMS procedures—Case 2</caption><thead><tr><th rowspan="3">Mode</th><th colspan="12">Number of constrained precessional modes retained</th></tr><tr><th colspan="3">20</th><th colspan="3">12</th><th colspan="3">8</th><th colspan="3">4</th><th colspan="3">2</th></tr><tr><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th><th>Real modes</th><th>Complex modes</th></tr></thead><tbody><tr><td>1B I</td><td>224.13</td><td>224.13</td><td>224.13</td><td>224.13</td><td>224.12</td><td>224.12</td><td>224.16</td><td>224.2</td><td>224.27</td><td>224.2</td><td>224.16</td><td>224.2</td><td>224.27</td><td>224.2</td></tr><tr><td>R</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-8.04</td><td>-7.987</td><td>-8.015</td><td>-8.04</td><td>-8.04</td><td>-7.987</td><td>-8.015</td></tr><tr><td>1F I</td><td>226.05</td><td>226.05</td><td>226.05</td><td>226.05</td><td>226.05</td><td>226.05</td><td>225.63</td><td>225.9</td><td>224.4</td><td>225.03</td><td>225.63</td><td>225.9</td><td>224.4</td><td>225.03</td></tr><tr><td>R</td><td>-45.87</td><td>-45.87</td><td>-45.9</td><td>-45.9</td><td>-45.87</td><td>-45.87</td><td>-45.8</td><td>-45.9</td><td>-45.7</td><td>-45.89</td><td>-45.87</td><td>-45.9</td><td>-45.7</td><td>-45.89</td></tr><tr><td>2B I</td><td>988.32</td><td>988.32</td><td>988.32</td><td>988.32</td><td>988.37</td><td>988.32</td><td>987.17</td><td>987.98</td><td>1052</td><td>1041.3</td><td>988.32</td><td>987.17</td><td>1052</td><td>1041.3</td></tr><tr><td>R</td><td>-277.95</td><td>-277.95</td><td>-277.97</td><td>-278.02</td><td>-277.91</td><td>-277.91</td><td>-269.2</td><td>-267.83</td><td>-241.5</td><td>-256.37</td><td>-277.95</td><td>-277.91</td><td>-269.2</td><td>-256.37</td></tr><tr><td>2F I</td><td>1031.0</td><td>1031.0</td><td>1031.03</td><td>1031.03</td><td>1031.0</td><td>1031.0</td><td>1031.0</td><td>1031.5</td><td>1087.1</td><td>1076.8</td><td>1031.0</td><td>1031.0</td><td>1087.1</td><td>1076.8</td></tr><tr><td>R</td><td>-205.31</td><td>-205.31</td><td>-205.35</td><td>-205.35</td><td>-205.3</td><td>-205.3</td><td>-199.9</td><td>-198.32</td><td>-180.7</td><td>-192.65</td><td>-205.31</td><td>-205.35</td><td>-199.9</td><td>-192.65</td></tr><tr><td>3B I</td><td>1893.95</td><td>1893.95</td><td>1894.4</td><td>1894.4</td><td>1895.3</td><td>1895.3</td><td>1894.3</td><td>1897.2</td><td>2347.7</td><td>2288.2</td><td>1893.95</td><td>1893.95</td><td>1894.4</td><td>1894.4</td></tr><tr><td>R</td><td>-44.3</td><td>-44.3</td><td>-44.23</td><td>-44.23</td><td>-44.3</td><td>-44.3</td><td>-43.7</td><td>-48.5</td><td>-117.4</td><td>-79.43</td><td>-44.3</td><td>-44.3</td><td>-48.5</td><td>-79.43</td></tr><tr><td>3F I</td><td>1922.7</td><td>1922.7</td><td>1923.2</td><td>1923.2</td><td>1923.4</td><td>1924.1</td><td>1923</td><td>1928.4</td><td>2411.5</td><td>2333.1</td><td>1922.7</td><td>1922.7</td><td>1923.2</td><td>1923.2</td></tr><tr><td>R</td><td>-45.77</td><td>-45.77</td><td>-45.7</td><td>-45.69</td><td>-45.69</td><td>-45.8</td><td>-45.2</td><td>-49.6</td><td>-117.7</td><td>-95.92</td><td>-45.77</td><td>-45.77</td><td>-45.69</td><td>-45.92</td></tr><tr><td>4B I</td><td>4233.1</td><td>4233.1</td><td>4233.4</td><td>4233.4</td><td>4233.4</td><td>4236.0</td><td>5253.4</td><td>5088</td><td>—</td><td>—</td><td>4233.1</td><td>4233.1</td><td>4233.4</td><td>4233.4</td></tr><tr><td>R</td><td>-158.48</td><td>-158.48</td><td>-158.48</td><td>-158.48</td><td>-158.48</td><td>-159.7</td><td>-76.7</td><td>-130.05</td><td>—</td><td>—</td><td>-158.48</td><td>-158.48</td><td>-159.7</td><td>-130.05</td></tr><tr><td>4F I</td><td>4305.32</td><td>4305.32</td><td>4305.52</td><td>4305.52</td><td>4305.52</td><td>4308.9</td><td>5345.2</td><td>5108.2</td><td>—</td><td>—</td><td>4305.32</td><td>4305.32</td><td>4305.52</td><td>4305.52</td></tr><tr><td>R</td><td>-153.89</td><td>-153.89</td><td>-153.85</td><td>-153.85</td><td>-155.2</td><td>-152.2</td><td>-73.0</td><td>-125.62</td><td>—</td><td>—</td><td>-153.89</td><td>-153.89</td><td>-153.85</td><td>-153.85</td></tr></tbody></table>

Notes: Real = real modes, Com = complex modes, I, R = imaginary and real parts of complex eigenvalues (rad/s).
<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic configuration of rotor bearing for Case 1—eigenvalue analysis using real and complex mode CMS.</figcaption></figure>

used if there is a necessity in terms of accuracy. In the following section, a real mode CMS procedure is used to calculate imbalance response, stability, and bifurcations of a journal bearing supported flexible rotor systems.

<figure><img src="image_3.png" /><figcaption>Figure 2. Schematic configuration of rotor bearing for Case 2—eigenvalue analysis using real and complex mode CMS.</figcaption></figure>
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. (a) Model of a flexible rotor supported on identical journal bearings used for CMS analysis—Case 3. (b) Journal bearing notation.</figcaption></figure>

## 6. UNBALANCE RESPONSE AND STABILITY OF FLEXIBLE ROTORS SUPPORTED ON JOURNAL BEARINGS

Most rotor systems have some mass imbalance and this causes an excitation at the running speed. Hence, it is important to be able to predict the unbalance response of a
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Unbalance response at the (a) bearing and (b) disk (node 3)—flexible rotor supported on journal bearings (Case 3): ●, 2 modes; —, 3 modes; ---, unstable.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 5. (a) Leading Floquet multiplier versus speed and (b) imaginary part versus real part of lead Floquet multiplier—flexible rotor supported on journal bearings (Case 3): ●, 2 modes; —, 3 modes.</figcaption></figure>

MDOF rotor bearing system for different speed and other parameter ranges of interest. The shooting and continuation scheme in combination with the real modes CMS procedure that has been developed for this study is used here to calculate the unbalance response of a MDOF rotor with nonlinear supports. To illustrate the power and usefulness of the combined numerical scheme, an example rotor system supported on plain journal bearings is considered (Figure 3). An unbalance mass of $3 \times 10^{-5} \text{ kg-m}$ is assumed to be present at nodes 3 and 4. The same rotor that was used before in Case 1 is analyzed here. No bearing other than the journal bearings are used to support the load. The journal
bearings are located at the end nodes 1 and 6. The journal bearing parameters used by Glasgow and Nelson [18] are used here and are length = 2.85 cm (1.12 in.), diameter = 5.72 cm (2.25 in.), viscosity = 6.9 cP, radial clearance = 0.051 mm (0.002 in.),

<figure><img src="image_9.png" /><figcaption>Figure 6. X-Y orbit at the (a) bearing and (b) disk (node 3) for S = 1.5—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7. X-Y orbit at the (a) bearing and (b) disk (node 3) for S = 1.7—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8. X-Y orbit at the (a) bearing and (b) disk (node 3) for S = 2.0—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9. X-Y orbit at the (a) bearing and (b) disk (node 3) for S = 2.5—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
static bearing load = 395 N (88.9 lb). A short bearing model is used to calculate the fluid-film forces. A short bearing model [20] is used to calculate the pressure field and fluid-film forces. The expression for the pressure distribution in this case is given by

$$P(\theta, t) = -\frac{1}{H^3}\left(\frac{dH}{dt} + 0.5 \frac{dH}{d\theta}\right), \qquad (32)$$

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

<figure><img src="image_15.png" /><figcaption>Figure 10. Frequency spectrum at the (a) bearing and (b) disk for S = 1.5—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 11. Frequency spectrum at the (a) bearing and (b) disk for S = 1.7—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
where fluid-film thickness $H = 1 + X(t) \cos(\theta) + Y(t) \sin(\theta)$, in which $\theta$ is measured counterclockwise from the negative X-axis [Figure 3(b)]. The bearing forces non-dimensionalized with respect to inertia force, $mc\omega^2$, can be written as

$$F_x = \frac{1}{S} \frac{\mu R}{m\omega_n} \left(\frac{L}{C}\right)^3 I_1, \quad F_y = \frac{1}{S} \frac{\mu R}{m\omega_n} \left(\frac{L}{C}\right)^3 I_2, \qquad (33)$$

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

<figure><img src="image_19.png" /><figcaption>Figure 12. Frequency spectrum at the (a) bearing and (b) disk for $S = 2.0$—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 13. Frequency spectrum at the (a) bearing and (b) disk for S = 2.5—flexible rotor supported on journal bearings (Case 3).</figcaption></figure>

$$I_1 = \int_0^{2\pi} P(\theta, t) \cos(\theta) \, d\theta, \quad I_2 = \int_0^{2\pi} P(\theta, t) \sin(\theta) \, d\theta. \qquad (34)$$

The equations of motion were non-dimensionalized as follows:

$$X = x/c, \quad Y = y/c, \quad X' = \dot{x}/(\omega c), \quad Y' = \dot{y}/(\omega c), \quad X'' = \ddot{x}/(c\omega^2), \\ Y'' = \ddot{y}/(c\omega^2), \quad \tau = \omega t.$$
The non-dimensionalization $\tau = \omega t$ makes it possible to integrate the equations of motion over a period of revolution $2p$ (as required by the shooting method) irrespective of the spin speed $\omega$. With this, the system can be quickly analyzed for various spin speeds. Also, the non-dimensionalization is useful for FFT calculations since $\omega$ does not directly appear in the calculations. The non-dimensional equations were utilized with the CMS and shooting/continuation scheme and the results for the journal bearing case are given below.

Figures 4(a) and (b) show the amplitude response of the rotor at the bearing and disk locations (nodes 1 and 3), calculated using two and three precessional modes in the combined CMS shooting procedure. The amplitude of the periodic response with the period the same as the running speed (period-1) is plotted over a wide, non-dimensional, speed range. The speed has been non-dimensionalized with respect to the nominal natural frequency $\omega_n$ of the rotor system calculated assuming the rotor is simply supported. The original rotor model has 24 degrees of freedom. Using the real mode CMS procedure, the number of degrees of freedom has been reduced to six (4 boundary co-ordinates + 2 precessional modes) and seven (4 boundary coordinates + 3 precessional modes) for the two modes and three modes cases, respectively. The figures show that the results obtained using two precessional modes is almost identical to those calculated using three modes. The dashed lines indicate unstable solution. The period-1 solution becomes unstable when the spin speed exceeds 1.65 times the nominal natural frequency ($S = 1.65$) and bifurcates into a period-2 solution. This is seen from the plot of the Floquet multipliers in Figures 5(a) and (b). Figure 5(a) shows the magnitude of the leading Floquet multiplier with respect to the non-dimensional speed. When the speed is 1.65, the magnitude of the leading Floquet multiplier is unity and at higher speeds this value is above unity indicating that the period-1 solution is unstable above $S = 1.65$. The type of bifurcation can be understood from Figure 5(b). It is seen that the Floquet multipliers exit the unit circle through the negative real axis indicating that the bifurcation is a period-doubling bifurcation. Such bifurcation results in the birth of a subharmonic vibration with a frequency that is exactly one-half the running speed. This subsynchronous vibration is called oil-whirl by rotordynamicists, based on linear analysis. Figures 6(a) and (b) show the orbital motion of the rotor at the bearing and unbalance location (nodes 1 and 3) before the period-doubling bifurcation (running speed less than the threshold speed). Both orbits show a periodic motion of period-1. The motion at the unbalance location is larger than that at the bearing. Figures 7(a) and (b) show the motion at nodes 1 and 3 at a non-dimensional speed of 1.7 which is just above the threshold speed 1.65. The orbits are much larger than the ones below the threshold speed and show a small convolution suggesting the presence of another frequency. (The loop typical of subsynchronous motion cannot be deciphered from the plot as drawn due to the presence of dominant subsynchronous frequency.) It is seen that the motion is approximately centered near the bearing center. As the speed increases, the orbit sizes also increase [Figures 8(a, b), 9(a, b)]. Figures 10(a) and (b) show the frequency spectra at nodes 1 and 3 for $S = 1.5$ (below the threshold speed). The spectra show a predominant response at the running speed which is due to unbalance. Figures 11(a) and (b) show the frequency spectrum at a speed parameter 1.7 which is just above the threshold speed. A dominant response at half-running speed is seen. This subsynchronous motion is due to the period-doubling bifurcation as seen from the Floquet multiplier plots before. Figures 12(a, b), 13(a, b) show the frequency spectra at speeds $S = 2.0$ and 2.5 respectively. These plots reveal that as the speed is increased the subsynchronous response tracks the running speed at exactly half that speed and increases in amplitude with speed.
# 7. CONCLUSIONS

The results of this study show that real component modes can be used to reduce the order of large order systems with local non-linearities without losing essential dynamics of the original system. The accuracy of modal parameters predicted using real modes CMS is comparable to those obtained using complex modes CMS (Tables 1 and 2). The real modes CMS procedure, in combination with a shooting/continuation scheme, was used to study the unbalance response and stability of a 24 DOF rotor supported on journal bearings and the results show that the system could be reduced to 6 DOF (4 boundary co-ordinates + the 2 lowest precessional modes) by retaining just 2 constrained normal modes, and the reduced system essentially displays all the dynamics of the original system. Truncation of higher modes did not have any significant effect on the threshold speed of instability. The results were verified by adding more modes to assure convergence. The method of real modes CMS with shooting/continuation shows promise in analyzing transient, steady state, stability, and bifurcational behavior of large order nonlinear systems with strong non-linearities.

# REFERENCES

1. S. YAMANCHI 1983 *JSME* **49**, 1862–1868. The nonlinear of flexible rotors, 1st report, development of a new analysis technique.

2. Y. S. CHOI and S. T. NOAH 1987 *ASME Journal of Vibration, Acoustics, Stress and Reliability in Design* **109**, 255–261. Nonlinear steady state response of a rotor-support system.

3. Y. B. KIM and S. T. NOAH 1991 *ASME Journal of Applied Mechanics* **58**, 545–553. Stability and bifurcation analysis of oscillators with piecewise-linear characteristics: a general approach.

4. Y. B. KIM and S. T. NOAH 1991 *Nonlinear Dynamics* **2**, 215–234. Response and bifurcation analysis of MDOF rotor system with a strong local nonlinearity.

5. T. N. SHIAU and A. N. JEAN 1989 12th *Biennial ASME Conference on Mechanical Vibration and Noise*. Prediction of steady state response of flexible rotor systems with nonlinear supports.

6. Y. S. CHOI and S. T. NOAH 1992 *Nonlinear Dynamics* **3**, 105–121. Response and stability analysis of piecewise-linear oscillators under multi-forcing frequencies.

7. L. J. MACNEAL and E. J. HAHN 1983 *ASME Journal of Lubrication Technology* **105**, 22–28.

8. C. NATARAJ and H. D. NELSON 1989 *ASME Journal of Vibration, Acoustics, Stress, and Reliability in Design* **111**, 187–193. Periodic solutions in rotordynamic systems with nonlinear supports.

9. A. N. JEAN and H. D. NELSON 1990 *Journal of Sound and Vibration* **143**, 473–489. Periodic response investigation of large nonlinear rotordynamic systems using collocation.

10. L. J. MACLEAN and E. J. HAHN 1983 *ASME Journal of Lubrication of lubrication Technology* **105**, 22–28. Unbalance behavior of squeeze-film damped multi-mass flexible rotor bearing systems.

11. R. H. B. FEY, D. H. VAN CAMPEN and A. DE KRAKER 1992 *Nonlinear Vibrations, Proceedings, ASME Symposium on Nonlinear Vibrations* **144**, 159–166. Long term structural dynamics of mechanical systems with local nonlinearities.

12. P. SUNDARARAJAN and S. NOAH 1994 5th *Conference on Nonlinear Vibrations, Stability and Dynamics of Structures*. Subcritical Hopf bifurcation in rotor supported on finite-length journal bearing.

13. P. SUNDARARAJAN, S. NOAH and L. A. SAN ANDRES 1994 *Proceedings, International Federation of the Theory of Machines and Mechanisms (IFTOMM) Fourth International Conference on Rotordynamics* 333–340. Fluid inertia effects on the nonlinear response of a squeeze-film supported rigid rotor system.

14. P. SUNDARARAJAN and S. T. NOAH 1995 *ASME Journal of Vibration and Acoustics*. Dynamics of forced nonlinear systems using shooting methods with arc length continuation—application to rotor systems.

15. P. SUNDARARAJAN and S. T. NOAH 1995 Analytical study of the nonlinear response of squeeze-film damper supported rotors and comparison with experimental results—eccentric operation.
16. H. B. KELLER 1987 *Lecture Notes on Numerical Methods in Bifurcation Problems*. Tata Institute of Fundamental Research. Bombay: Springer.

17. S. CRANDALL and N. YEH 1987 *Rotating Machinery Dynamics, Proceedings, 11th Vibration Conference 1*, 79-84. Automtic generation of component modes for rotordynamic structures.

18. D. A. GLASGOW and H. D. NELSON 1980 *Transactions of ASME, Journal of Mechanical Design* **102**, 352-358. Stability analysis of rotor-bearing systems using component mode synthesis.

19. H. D. NELSON, S. J. RUSSEL and W. L. MEACHAM 1981 Transient response of rotor-bearing systems using component mode synthesis, part I—mathematical development.

20. D. CHILDS 1993 *Turbomachinery Rotordynamics—Phenomena, Modeling and Analysis*. New York: Wiley.

21. HOU, SHOU-NIEN 1969 *Shock and Vibration Bulletin* **40**, 25-39. Review of modal synthesis techniques and a new approach.

22. R. R. CRAIG JR. and Y. CHUNG 1981 *AIAA/ASME 20th Structures, Structural Dynamics and Material Conference, AIAA* 81-0560, 254-266. A generalized substructure coupling procedure for damped systems.

23. R. R. CRAIG JR. 1987 *Journal of Analytical and Experimental Modal Analysis* **2**, 59-72. A review of time-domain and frequency-domain component-mode synthesis methods.

24. W. C. HURTY 1965 *AIAA Journal* **3**, 678-685. Dynamic analysis of structural systems using component modes.

25. R. R. CRAIG JR. and M. C. BAMPTON 1968 *AIAA Journal* **6**, 1313-1319. Coupling of substructures for dynamic analysis.

26. D. F. LI and E. J. GUNTER 1982 *ASME Journal of Engineering for Power* **104**, 552-560. Component mode synthesis of large rotor systems.

27. M. L. ADAMS 1980 *Journal of Sound and Vibration* **71**, 129-144. Nonlinear dynamics of flexible multi-bearing rotors.

28. R. H. MACNEAL 1971 *Journal of Computers and Structures* **1**, 581-601. A hybrid method of component mode synthesis.

29. S. RUBIN 1975 *AIAA Journal* **13**, 995-1006. An improved component-mode representation.

30. R. SEYDEL 1988 *From Equilibrium to Chaos, Practical Bifurcation and Stability Analysis*. Amsterdam: Elsevier Science.

31. T. S. PARKER and L. O. CHUA 1989 *Practical Numerical Algorithms for Chaotics Systems*. New York: Springer.

32. BALACHANDRAN and A. NAYFEH 1995 *Applied Nonlinear Dynamics*. New York: Wiley.

33. M. ALUKO and H. CHANG 1984 *Computers and Chemical Engineering* **8**(6), 355-365. PEFLOQ; an algorithm for the bifurcational analysis of periodic solutions of autonomous systems.

34. P. SUNDARARAJAN and S. T. NOAH 1996 *ASME Journal of Vibration and Acoustics*. Dynamics of forced nonlinear systems using shooting methods with arc length continuation—application to rotor systems.

35. E. DOEDEL, H. B. KELLER and J. P. KERNEVEZ 1991 *International Journal of Bifurcation and Chaos* **1**, 745-752.

36. M. KUBICEK and V. HLAVCEK 1983 *Numerical Solution of Nonlinear Boundary Value Problems with Applications*. Englewood Cliffs, N.J.: Prentice-Hall.

## NOMENCLATURE

*   *A* real eigenvector matrix of the system

*   *c* radial clearance in journal bearing/squeeze film damper = *R* − *r*

*   CMS Component Mode Synthesis

*   *C* component damping matrix

*   *E* α<sup>*T*</sup>A<sub>1</sub>α

*   *e* eccentricity

*   *F* α<sup>*T*</sup>B<sub>1</sub>α

*   *f* static load per bearing

*   *f*<sub>*x*</sub>, *f*<sub>*y*</sub> fluid forces in inertial *X* and *Y* frames

*   *F*<sub>*x*</sub>, *F*<sub>*y*</sub> non-dimensional fluid forces in inertial *X* and *Y* frames = *f*<sub>*x*</sub>/mcω<sup>2</sup>, *f*<sub>*y*</sub>/mcω<sup>2</sup>

*   *h* ( *ẋ*<sub>*b*</sub>, *q̇*, *x*<sub>*b*</sub>, *q*)<sup>*T*</sup>

*   *I* unit matrix

*   *L* bearing length
$M$ component mass matrix

$K$ component stiffness matrix

$P(\theta, t)$ pressure distribution in the damper

$Q$ $\alpha^T F_1 \alpha$

$q$ modal coordinates in CMS procedure

$r, R$ radii of the journal and bearing respectively

$S$ speed parameter = $\omega/\omega_n$ (defined for journal bearing systems)

$t$ tme, $\tau = \omega t$.

$u$ imbalance; $U = u/c$

$x, y$ co-ordinates of journal center with respect to X-Y system fixed to bearing center

$X$ $x/c$

$Y$ $y/c$

$y$ $(\dot{x}_b, x_b, \dot{q}, q)^T$

$\epsilon$ $e/c$

$\alpha$ transformation matrix from vector $y$ to vector $h$

$\omega$ angular spin speed of the rotor

$\omega_n$ natural frequency = $\sqrt{k/m}$

$\mu$ dynamic viscosity of bearing lubricant

$\theta$ angular position measured with respect to negative X axis (short bearings)

$\sigma$ somerfeld number = $\mu\omega LR^3/fc^2$

## Subscript

$0$ static equilibrium value of a variable (for journal bearings)

$i$ internal degrees of freedom in CMS procedure

$b$ boundary degrees of freedom in CMS procedure

$nsb$ non-symmetric, boundary DOF

$nsi$ non-symmetric, internal DOF

$s$ symmetric