# Nonlinear Normal Modes of a Rotating Shaft Based on the Invariant Manifold Method

Mathias Legrand, Dongying Jiang, Christophe Pierre, Steven Shaw

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

► To cite this version:

Mathias Legrand, Dongying Jiang, Christophe Pierre, Steven Shaw. Nonlinear Normal Modes of a Rotating Shaft Based on the Invariant Manifold Method. International Journal of Rotating Machinery, 2002, 10.1155/S1023621X04000338. hal-01315570

HAL Id: hal-01315570

https://hal.science/hal-01315570

Submitted on 13 May 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear Normal Modes of a Rotating Shaft Based on the Invariant Manifold Method

**Mathias Legrand**

*GeM, Equipe Structures et Simulations, Ecole Centrale de Nantes, Nantes, France*

**Dongying Jiang**

*Department of Mechanical Engineering, University of Michigan, Ann Arbor, Michigan, USA*

**Christophe Pierre**

*Department of Mechanical Engineering, University of Michigan, Ann Arbor, Michigan, USA*

**Steven W. Shaw**

*Department of Mechanical Engineering, Michigan State University, East Lansing, Michigan, USA*

The nonlinear normal mode methodology is generalized to the study of a rotating shaft supported by two short journal bearings. For rotating shafts, nonlinearities are generated by forces arising from the supporting hydraulic bearings. In this study, the rotating shaft is represented by a linear beam, while a simplified bearing model is employed so that the nonlinear supporting forces can be expressed analytically. The equations of motion of the coupled shaft-bearings system are constructed using the Craig-Bampton method of component mode synthesis, producing a model with as few as six degrees of freedom (d.o.f.). Using an invariant manifold approach, the individual nonlinear normal modes of the shaft-bearings system are then constructed, yielding a single-d.o.f. reduced-order model for each nonlinear mode. This requires a generalized formulation for the manifolds, since the system features damping as well as gyroscopic and nonconservative circulatory terms. The nonlinear modes are calculated numerically using a nonlinear Galerkin method that is able to capture large amplitude motions. The shaft response from the nonlinear mode model is shown to match extremely well the simulations from the reference Craig-Bampton model.

**Keywords** Invariant manifold method, Nonlinear normal modes, Oil-whip, Oil-whirl, Shaft-bearing system

Many rotating structures are supported by devices that are inherently nonlinear, e.g., journal bearings. The dynamic analysis of nonlinear rotating systems has been the subject of a number of studies, for example, those of Yamauchi (1983) or Kim and Noah¹ (1991a), who employed the method of harmonic balance. Choi and Noah (1987) added discrete Fourier transform procedures to the harmonic balance method and also included subharmonic response components. Kim and Noah² (1991b) used dynamic condensation techniques in conjunction with the harmonic balance method, in order to reduce the size of the system models. The objective of the present work is to develop reduced-order models of nonlinear shaft-bearings systems using the invariant manifold-based nonlinear normal mode methodology.

Nonlinear normal modes (NNM) provide a general framework for the construction of reduced-order models for nonlinear systems. The concept of NNM was first introduced by Rosenberg (1966) with the study of conservative, symmetric, nonlinear systems. A NNM was defined in the configuration space so that its application was strictly limited to systems without gyroscopic effects and damping. Shaw and Pierre (1993) extended the definition of nonlinear normal modes using invariant manifold techniques, wherein a NNM is defined as a two-dimensional invariant surface in the phase space, which is tangent to the hyperplane that represents the corresponding mode of the linearized model. The NNM response is thus captured by a single-d.o.f nonlinear oscillator. Using this more general definition, a systematic construction method for NNM has been proposed by Boivin (1995) and Pescheck (2000) for nonlinear systems with quadratic and cubic nonlinearities, including systems with a large number of d.o.f. and for large amplitude motions. Nayfeh and Nayfeh
(1994) also used the invariant manifold approach to construct NNM using perturbation methods for weakly nonlinear systems.

In this article, the NNM are constructed for an established model of a nonlinear shaft-bearings system, which consists of a linear rotating shaft supported by short nonlinear journal bearings at its two ends. Since the linearized system is gyroscopic, damped, and (nonconservative) circulatory—thus featuring a non-symmetric stiffness matrix—an extension of the invariant manifold approach is required in order to accommodate these effects. The NNM are shown to provide very accurate reduced-order models of the shaft-bearings system.

The article is organized as follows. In the first section, the mathematical model of the shaft-bearings system is derived using component mode synthesis (CMS). Then, in the second section, the vibration modes of the linearized system are investigated. Finally, in the third section, the theory of NNM is extended to systems whose linearized counterpart is gyroscopic, nonconservative, and circulatory, and the individual NNM invariant manifolds of the nonlinear shaft-bearings system are calculated.

## THE ROTATING SHAFT-BEARINGS SYSTEM

A diagram of the system of interest, a rotating shaft supported by short journal bearings at its two ends, is shown in Figure 1. A Rayleigh beam with uniform cross-section properties is used to model the shaft, which is assumed to be perfectly balanced. The forces created by the oil film in the bearings are nonlinear and can be represented as nonlinear boundary conditions for the beam.

In Figure 1, the inertial frame, $R_{XYZ}$, is fixed in space and the Y-axis passes through the centers of the bearings at the two ends. The shaft is defined by its length $L = 1$ m, outer diameter $D_2 = 0.0592$ m, and inner diameter $D_1 = 0.02$ m. The nominal clearance between the bearing and the shaft is set as $c = (D_3 - D_2)/2 = 5.1 \cdot 10^{-5}$ m, where $D_3$ is the inner diameter for both bearings. The length of the bearing is $L_b = 0.0285$ m, and the dynamic viscosity of the oil film is chosen as $\mu = 0.0068 \text{ N} \cdot \text{s/m}^2$.

<figure><img src="image_3.png" /><figcaption>FIGURE 1
Schematic of the shaft-bearings system.</figcaption></figure>

## The Shaft Model

From the Rayleigh beam theory, the kinetic energy, $T$, and strain energy, $U$, of the rotating shaft are given by, respectively:

$$ T = \frac{\rho A}{2} \int_{0}^{L} (\dot{u}^2 + \dot{w}^2) dy + \frac{\rho I}{2} \int_{0}^{L} \left( \left( \frac{\partial \dot{u}}{\partial y} \right)^2 + \left( \frac{\partial \dot{w}}{\partial y} \right)^2 \right) dy \\ + \rho I \Omega^2 - 2\rho I \Omega \int_{0}^{L} \frac{\partial \dot{u}}{\partial y} \frac{\partial \dot{w}}{\partial y} dy \quad [1] $$

$$ U = \frac{EI}{2} \int_{0}^{L} \left( \left( \frac{\partial^2 u}{\partial y^2} \right)^2 + \left( \frac{\partial^2 w}{\partial y^2} \right)^2 \right) dy \quad [2] $$

where $u(y, t)$ and $w(y, t)$ are the displacements of the points on the neutral axis of the shaft in the X and Z directions, respectively, $\Omega$ is the constant angular velocity of the shaft, $A$ is the shaft's cross-sectional area, and $I$ is the second area moment of inertia of the cross section. The material parameters for the shaft are its Young's modulus $E = 2.1 \cdot 10^{11} \text{ Pa}$ and mass density $\rho = 7800 \text{ kg/m}^3$.

In order to obtain an efficient discretized model for the rotating shaft, the method of component mode synthesis developed by Craig and Bampton (1968) is applied, where the shaft is the linear substructure constrained at its ends by the nonlinear bearings. This allows for the nonlinear effects of the supporting bearings to be captured solely by the d.o.f. corresponding to the static constraint modes. The shaft displacements $u$ and $w$ are thus each expanded as a linear combination of the modes of free vibration of the shaft pinned at both ends and two constraint modes, each corresponding to a rigid body motion of the shaft induced by a unit displacement at one of its ends:

$$ \begin{aligned} u(y, t) &= \sum_{i=1}^{m} \Phi_i(y) a_i(t) \\ w(y, t) &= \sum_{i=1}^{m} \Phi_i(y) b_i(t) \end{aligned} \quad [3] $$

where

$$ \Phi_i(y) = \begin{cases} \sin\left(\frac{i\pi y}{L}\right) & i = 1, \dots, m-2 \\ \frac{y}{L} & i = m-1 \\ 1 - \frac{y}{L} & i = m \end{cases} $$

The first $(m-2)$ expansion functions are the mode shapes of the simply supported shaft and the final two are the static constraint modes. We substitute the expansion functions for $u$ and $w$ into Hamilton's principle,

$$ \int_{t_1}^{t_2} (\delta T - \delta U + \delta W) dt = 0 \quad \forall (t_1, t_2) \quad [4] $$

where $\delta W$ is the virtual work done by the gravity force and the nonlinear supporting forces of the two bearings. The resulting discretized model is given as:

$$ [M_1]\{\ddot{x}\} + [G_1]\{\dot{x}\} + [K_1]\{x\} = \{F\} \quad [5] $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2
Schematic of the journal bearing model.</figcaption></figure>

where

$$ \{x\} = \begin{Bmatrix} \{a\} \\ \{b\} \end{Bmatrix}, [M_1] = \begin{pmatrix} \rho S \begin{bmatrix} [A] & [0] \\ [0] & [A] \end{bmatrix} + \rho I \begin{bmatrix} [B] & [0] \\ [0] & [B] \end{bmatrix} \end{pmatrix}, \\ [G_1] = 2\rho I \Omega \begin{bmatrix} [0] & [B] \\ -[B] & [0] \end{bmatrix}, [K_1] = EI \begin{bmatrix} [C] & [0] \\ [0] & [C] \end{bmatrix} $$

Here {$a$} and {$b$} are m-vectors of the generalized coordinates generated by the CMS formulation, $a_i$ and $b_i$, respectively. The matrices [A], [B], and [C] are listed as follows:

$$ \begin{align*} [A] &= \int_{0}^{L} [\Phi(y)]^{T} [\Phi(y)] dy \\ [B] &= \int_{0}^{L} [\Phi_{,y}(y)]^{T} [\Phi_{,y}(y)] dy \\ [C] &= \int_{0}^{L} [\Phi_{,yy}(y)]^{T} [\Phi_{,yy}(y)] dy \end{align*} $$

The matrices {$[M_1]$, $[G_1]$, and $[K_1]$}, are the inertia, gyroscopic, and stiffness matrices, respectively. Note that {$x$} is defined with respect to the inertial frame $R_{XYZ}$. The force vector {$F$} has the form:

$$ \{F\} = -\rho g A \int_{0}^{L} \{\Phi(y)\} dy + \{F_{2}(x, \dot{x})\} = \{F_{1}\} + \{F_{2}(x, \dot{x})\} \quad [6] $$

## The Bearing Forces

In Figure 2, $R_{XYZ}$ is the inertial frame and $\Omega$ is the angular velocity, defined in Figure 1. The nominal thickness of the oil film is $h$, $e$ is the eccentricity between the bearing axis and the shaft axis, and $\phi$ is the attitude angle of the line connecting the bearing and shaft centers with respect to the Z-axis. The horizontal and vertical displacements of the center of the shaft (the journal) in the bearing are denoted as $x_j$ and $z_j$, respectively (j for journal). The thickness of the oil film, $h$, can be expressed as:

$$ h = c - z_j \cos(\theta + \phi) + x_j \sin(\theta + \phi) \quad [7] $$

where $c$ is the nominal clearance between the shaft and the bearing. Based on Reynolds' equation, the pressure of the oil film can be modeled as:

$$ \frac{\partial}{\partial y} \left( \frac{h^3}{6\mu} \frac{\partial p}{\partial y} \right) + \frac{1}{R^2} \frac{\partial}{\partial \theta} \left( \frac{h^3}{6\mu} \frac{\partial p}{\partial \theta} \right) = \Omega \frac{\partial h}{\partial \theta} + 2 \frac{\partial h}{\partial t} \quad [8] $$

In Equation (8), $\mu$ is the fluid viscosity, $R$ is the outer radius of the beam, and $p$ is the fluid film pressure. The short bearing assumption implies that $R^2$ is preponderant in Equation (8), so that the second term on the left-hand side in Reynolds' equation can be neglected. This yields:

$$ \frac{h^3}{6\mu} \frac{\partial^2 p}{\partial y^2} = \Omega(z_j \sin(\theta + \phi) + x_j \cos(\theta + \phi)) - 2(\dot{z}_j \cos(\theta + \phi) - \dot{x}_j \sin(\theta + \phi)) \quad [9] $$

The boundary conditions over y are $p(\theta, 0) = p(\theta, L_b) = 0$, where $L_b$ is the length of the bearing and $\theta$ is integrated over $[0, \pi]$ instead of $[0, 2\pi]$ due to cavitation effects described by Vance (1988). Then the two forces created by the pressure field are the integrals of the pressure over the fluid-film surface contact. A dimensionless analysis of the problem leads to the following definitions:

$$ \begin{align} Z_j &= \frac{z_j}{c}, & X_j &= \frac{x_j}{c}, & \dot{Z}_j &= \frac{\dot{z}_j}{\Omega c}, & \dot{X}_j &= \frac{\dot{x}_j}{\Omega c} \tag{10} \\ H &= 1 - Z_j \cos(\theta + \phi) + X_j \sin(\theta + \phi) & & \end{align} $$

where $H$ is the dimensionless fluid-film thickness. Then, the resultant forces $F_X$ and $F_Z$ in the X and Z directions can be obtained as (see Lee, 1993):

$$ F_X = -\frac{\mu R L^3 \Omega}{2c^2} \int_0^\pi \left( \frac{Z_j \sin(\theta + \phi) + X_j \cos(\theta + \phi) - 2(\dot{Z}_j \cos(\theta + \phi) - \dot{X}_j \sin(\theta + \phi))}{H^3} \right) \sin(\theta + \phi) d\theta \quad [11a] $$

$$ F_Z = \frac{\mu R L^3 \Omega}{2c^2} \int_0^\pi \left( \frac{Z_j \sin(\theta + \phi) + X_j \cos(\theta + \phi) - 2(\dot{Z}_j \cos(\theta + \phi) - \dot{X}_j \sin(\theta + \phi))}{H^3} \right) \cos(\theta + \phi) d\theta \quad [11b] $$

where the first term is the gravitational load, and the second term {$F_2$} is the supporting force at the journal bearing, which is considered next.

Both forces are nonlinear in the displacements and velocities of the shaft's ends, involving complicated integrals that can be
obtained analytically using commercial mathematical software such as Maple®.

## THE LINEARIZED MODEL

Complex modal analysis is applied to the system linearized about its equilibrium position to determine the natural modes of vibration and the shaft's first critical speed.

### The Equilibrium Position

The dimensionless equilibrium position of the shaft center is given by Ocvirck (1952) as:

$$ \begin{cases} Z_e = -\epsilon_e \cos \phi_e \\ X_e = \epsilon_e \sin \phi_e \end{cases} \quad \text{where } \phi_e = \arctan \left( \frac{\pi \sqrt{1 - \epsilon_e^2}}{4\epsilon_e} \right) \quad [12] $$

and the dimensionless eccentricity $\epsilon_e = e_e/c$ is the solution of:

$$ \frac{\mu N}{f} \frac{R}{2c^2} L_b^3 = \frac{(1 - \epsilon^2)^2}{\pi \epsilon \sqrt{\pi^2(1 - \epsilon^2) + 16\epsilon^2}} \quad [13] $$

where N (with $\Omega = 2\pi N$) is the angular frequency and f is the load, which in this study is given by the weight of the beam. [The locus of the equilibrium position versus $\Omega$ is shown on Fig. 3]. By linearizing the bearing forces with respect to this equilibrium

position, the stiffness matrix ($K_2$) (which is not symmetric) and the damping matrix ($D_2$) can be obtained as:

$$ K_{2IJ} = \left. \frac{\partial F_I}{\partial J} \right|_{\text{with } I = (X, Z) \text{ and } J = (X, Z)} = \frac{\partial F_I}{\partial J} $$

where $Z = Z_j - Z_e$ and $X = X_j - X_e$. The eight dimensionless coefficients for these matrices are given by Lee (1993) as follows:

$$ K_{ZZ} = - \frac{4(\pi^2 + (32 + \pi^2)\epsilon^2 + 2(16 - \pi^2)\epsilon^4)}{1 - \epsilon^2} \cdot Q(\epsilon) $$

$$ K_{ZX} = \frac{\pi(\pi^2 + (32 + \pi^2)\epsilon^2 + 2(16 - \pi^2)\epsilon^4)}{\epsilon\sqrt{1 - \epsilon^2}} \cdot Q(\epsilon) $$

$$ K_{XZ} = - \frac{\pi(\pi^2 - 2\pi^2\epsilon^2 - (16 - \pi^2)\epsilon^4)}{\epsilon\sqrt{1 - \epsilon^2}} \cdot Q(\epsilon) $$

$$ K_{XX} = -4(2\pi^2 + (16 - \pi^2)\epsilon^2) \cdot Q(\epsilon) $$

$$ C_{ZZ} = - \frac{2\pi(\pi^2 + 2(24 - \pi^2)\epsilon^2 + \pi^2\epsilon^4)}{\epsilon\sqrt{1 - \epsilon^2}} \cdot Q(\epsilon) $$

$$ C_{ZX} = 8(\pi^2 - 2(8 - \pi^2)\epsilon^2) \cdot Q(\epsilon) $$

$$ C_{XZ} = C_{ZX} $$

<figure><img src="image_5.png" /><figcaption>FIGURE 3
Locus of the equilibrium position versus $\Omega$.</figcaption></figure>
$$C_{XX} = - \frac{2\pi\sqrt{1 - \varepsilon^2}(\pi^2 - 2(8 - \pi^2)\varepsilon^2)}{\varepsilon} \cdot Q(\varepsilon)$$

$$Q(\varepsilon) = \frac{1}{(\pi^2(1 - \varepsilon^2) + 16\varepsilon^2)^{3/2}}$$

The relationships between the dimensionless and the dimensional coefficients are:

$$K_{IJ} = \frac{k_{ij} \cdot c}{f} \quad \text{and} \quad C_{IJ} = \frac{c_{ij} \cdot c \cdot \Omega}{f}$$

## The Linear Modes

In order to study the modes of free vibration of this system, Equation (5) is linearized with respect to the equilibrium position. The linear system can be written in the following general form:

$$[M]\{\ddot{x}\} + [C]\{\dot{x}\} + [K]\{x\} = \{0\} \quad [14]$$

where $[M] = [M_1]$, $[C] = [G_1] - [D_2]$, $[K] = [K_1] - [K_2]$, and $\{x\}$ is defined relative to the equilibrium position.

Complex modal analysis can be applied to the state space formulation of the system. Defining the velocity vector as $\{y\} = \{\dot{x}\}$, Equation (14) can be written in first-order form as:

$$\{\dot{z}\} = [D]\{z\} \quad [15]$$

where

$$\{z\} = \begin{Bmatrix} \{x\} \\ \{y\} \end{Bmatrix} \quad \text{and} \quad [D] = \begin{bmatrix} [O] & [I] \\ -[M]^{-1}[K] & -[M]^{-1}[C] \end{bmatrix}$$

Since the matrix [D] is real, those eigensolutions that are complex must occur in conjugate pairs. Since the NNM are to be defined in terms of real modal coordinates, a traditional complex mode transformation matrix would lead to difficulties in constructing the corresponding invariant manifolds. Thus, herein a real modal transformation is chosen that renders the matrix [D] block diagonal. More details on this procedure are given by Hirsch and Smale (1974). This is achieved by constructing the matrix $[Z']$, whose columns are the real eigenvectors and the real and the imaginary parts of the complex eigenvectors, as follows:

$$[Z'] = [\mathfrak{R}z_1, \mathfrak{R}z_2, \dots, \mathfrak{R}z_p, \mathfrak{R}z_p, z_{2p+1}, \dots, z_{2N}] \quad [16]$$

where N is the total number of d.o.f. for the rotating shaft system described by Equation (5), and is equal to 2 m.

Letting $\{z(t)\} = [Z']\{\eta(t)\}$, Equation (15) is transformed into:

$$\{\dot{\eta}(t)\} = [\Lambda']\{\eta(t)\}, \quad [17]$$

with

$$[\Lambda'] = \begin{bmatrix} \mathfrak{R}\lambda_1 & \mathfrak{R}\lambda_2 & & & & & & \\ -\mathfrak{R}\lambda_1 & \mathfrak{R}\lambda_2 & & & & & & \\ & \ddots & \mathfrak{R}\lambda_p & \mathfrak{R}\lambda_{p-1} & & & & \\ & & -\mathfrak{R}\lambda_p & \mathfrak{R}\lambda_{p-1} & \lambda_{2p+1} & & & \\ & & & 0 & \ddots & \mathfrak{R}\lambda_1 & \mathfrak{R}\lambda_2 & \\ & & & & & \lambda_{2N} & & \end{bmatrix}$$

The resulting first-order equations are at worst pairwise coupled. Those equations that are uncoupled are associated with the real eigenvalues and correspond to overdamped modes, whose dynamics occur in a one-dimensional linear subspace. Motions in these modes are decaying and non-oscillatory. The modes associated with the complex conjugate pairs of eigenvalues correspond to the 2 × 2 diagonal blocks of $[\Lambda']$ and are for underdamped, oscillatory modes. The dynamics of these modes occur in two-dimensional subspaces. Motions in the underdamped modes consist of decaying oscillations, which are described by pairs of first-order differential equations:

$$\begin{cases} \dot{s}_i = -\zeta_i \omega_i s_i + \omega_i \sqrt{1 - \zeta_i^2} t_i \\ \dot{t}_i = -\omega_i \sqrt{1 - \zeta_i^2} s_i - \zeta_i \omega_i t_i \end{cases} \quad i = 1, \dots, p \quad [18]$$

where $\omega_i$ is the undamped natural frequency of the mode and $\zeta_i$ its modal damping ratio, both of which are defined from $\mathfrak{R}\lambda_i = -\zeta_i\omega_i$ and $\mathfrak{R}\lambda_i = \omega_i\sqrt{1 - \zeta_i^2}$. Also, $s_i = \eta_{2i-1}$ and $t_i = \eta_{2i}$ are the components of $\{\eta\}$ associated with the ith underdamped mode.

## Stability of The Linear Model

In the Lyapunov sense, the system is unstable if the real part of any eigenvalue is positive. The linear motion associated with such an instability becomes infinitely large with time and thus linear theory breaks down. The eight coefficients of the linear model given by Lee (1993) are strongly linked to the stability of the system, especially the cross-coupled stiffness terms. In this study we have examined the stability of the shaft-bearings system using two, four, and six normal modes of vibration to describe the shaft motion, corresponding, respectively, to shaft-bearings models with six, eight, and ten d.o.f.s. The results for these three models have been found to be in close agreement, indicating that a single normal mode in each direction is sufficient to represent the shaft's motion. This makes sense since
<figure><img src="image_6.png" /><figcaption>FIGURE 4
Stability of the linear model: largest real part of the eigenvalues versus Ω.</figcaption></figure>

the shaft is balanced and thus the first beam mode is dominant. Therefore, the six-d.o.f. model is employed in the remainder of the paper. Figure 4 shows that the linearized system is unstable when the angular velocity of the shaft is larger than 375 rad/s.

## Linear Mode Motions

As the convergence study indicated, a six-d.o.f. model is sufficient to capture accurately the lower oscillating modes of the system. Recall that this model includes one free vibration normal mode and two static constraint modes for the shaft in each of the X and Z directions. The eigenanalysis was carried out in the state space in terms of the shaft angular velocity.

It was found that up to approximately Ω = 550 rad/s, the shaft-bearings system possesses four, one-dimensional overdamped modes, corresponding to real eigenvalues, and four, two-dimensional oscillatory modes, corresponding to complex eigenvalues. Two of the four overdamped modes become a single two-dimensional oscillatory mode as Ω grows beyond 550 rad/s, and the last two overdamped modes turn into a sixth oscillatory mode for Ω > 700 rad/s as shown in Figures 5a and 5b. In the following, we focus on the oscillatory motions at the rotation speed Ω = 100 rad/s. Since both real and imaginary parts are very close together for the two first modes, their undamped natural frequency and modal damping ratio are very close too, as depicted in Figures 6 and 7.

The physical nature of the modes can be depicted as follows. In Figures 6–9, curves are drawn for each oscillating mode in the X and Z coordinates that display the displacements of both ends of the beam as they move in the bearings during a modal motion. The accompanying sketches of the shaft show it as straight, whereas in reality it deflects in its first vibration mode (although its amplitude is relatively small compared to the bearing deflections).

## Comparison Between Nonlinear and Linear Models

Using a very small disturbance from the equilibrium position as initial conditions for both the linear and nonlinear models, the time simulation results are seen to be very close in Figure 10. For Ω = 100 rad/s, which lies in the stable range, after a short time the beam settles back to its equilibrium position.

As seen in Figure 11, for a much larger shaft rotation speed, Ω = 2550 rad/s, the linear model is unstable, as the angular velocity is then greater than the critical speed, Ω_c ≈ 375 rad/s, whereas the nonlinear model reaches a limit cycle corresponding to oil whirl detailed by Muszynska (1986). For a lightly loaded shaft, the whirl occurs at a frequency that is exactly half the shaft’s angular velocity, and here the whirl motion frequency is 1225 rad/s. The thick black line in Figure 11 represents the largest displacement of the center point of the shaft at the bearing location that is kinematically allowed. As shown in Figure 11, the oil whirl phenomenon cannot be predicted by linearized stability theory.
<figure><img src="image_7.png" /><figcaption>FIGURE 5a
The largest four real parts of the eigenvalues versus Ω.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 5b
Imaginary part of the eigenvalues versus Ω.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 6
First linear mode shape at $\Omega = 100$ rad/s: $\omega_1 = 50.647$, $\zeta_1 = 0.1166$.</figcaption></figure>

## NONLINEAR NORMAL MODES

For a linear system, an important property of a modal motion is that if one knows the motion of a single generalized coordinate, then the motions of all others coordinates are specified by that mode’s eigenvector. The same property is used for the nonlinear case, wherein a motion in a nonlinear normal mode is defined such that all displacements and velocities are functionally related to a single displacement-velocity pair. In a NNM the nonlinear system thus behaves essentially as a single-d.o.f. oscillator. Even though NNM are defined using a linear normal mode concept, the similarities stop there. Many results obtained by NNM have no counterparts in linear systems, such as amplitude-dependant frequencies and internal resonances (see Nayfeh and Mook [1979]), where two NNM exchange energy during the motion.

Using the concept of invariant manifolds, Shaw and Pierre (1993) defined a NNM as a motion that lies on a two-dimensional invariant manifold in the system’s phase space. Here “invariant” indicates that any motion initiated on the manifold will remain on it for all times. A single displacement-velocity pair is chosen as master coordinates, which characterize the individual nonlinear mode motion. All the remaining “slave” coordinates are parameterized by these two master coordinates and give the constrained conditions. Previous work on the invariant manifold method has been restricted to systems with diagonal stiffness and damping matrices in the modal coordinates space. However, in order to construct NNM for the rotating shaft system, a new formulation had to be introduced for obtaining the PDE governing the invariant manifold. This new formulation is more general and can be applied to systems with nonproportional damping
<figure><img src="image_10.png" /><figcaption>FIGURE 7
Second linear mode shape at $\Omega = 100$ rad/s: $\omega_2 = 50.676$, $\zeta_2 = 0.1112$.</figcaption></figure>

forces, gyroscopic effects, and nonsymmetric circulatory stiffness matrices.

## NNM Formulation

In this section we use Equation (5) written with respect to the equilibrium position. The equations of motion become:

$$[M]\{\ddot{x}\} + [C]\{\dot{x}\} + [K]\{x\} = \{F_{NL}(x, \dot{x})\} \quad [19]$$

and

$$\{F_{NL}(x, \dot{x})\} = \{F_2(x + x_e, \dot{x})\} - [K_2]\{x\} - [D_2]\{\dot{x}\} + \{F_1\} - [K_1]\{x_e\}$$

where $\{x\}$ is defined with respect to the equilibrium position $\{x_e\}$. All the matrices and forces in this equation have been defined in previous sections. Equation (19) are transformed following

the procedure of Equation (15) and the linear transformation of Equation (17) is carried out. It is assumed that the real linear transformation can be written as follows:

$$[\Lambda'] = \begin{bmatrix} \alpha_1 & \beta_1 \\ \beta_2 & \alpha_2 \\ \vdots & \vdots \\ \alpha_{2N-1} & \beta_{2N-1} \\ \beta_{2N} & \alpha_{2N} \end{bmatrix} \quad [20] \\ \begin{cases} \alpha_k = \alpha_{k+1} \\ \beta_k = -\beta_{k+1} \end{cases} \quad k = 1, 3, \dots, p-1, \\ \begin{cases} \alpha_k \neq \alpha_{k+1} \\ \beta_k = \beta_{k+1} = 0 \end{cases} \quad k = p+1, p+3, \dots, 2N-1$$
<figure><img src="image_11.png" /><figcaption>FIGURE 8
Third linear mode shape at $\Omega = 100$ rad/s: $\omega_3 = 775.495$, $\zeta_3 = 0.0457$.</figcaption></figure>

Then, the new equations of motion are:

$$ \begin{cases} \dot{s}_1 \\ \dot{t}_1 \\ \vdots \\ \dot{s}_N \\ \dot{t}_N \end{cases} = \begin{bmatrix} \alpha_1 & \beta_1 \\ \beta_2 & \alpha_2 \\ \vdots & \vdots \\ \alpha_{2N-1} & \beta_{2N-1} \\ \beta_{2N} & \alpha_{2N} \end{bmatrix} \cdot \begin{bmatrix} s_1 \\ t_1 \\ \vdots \\ s_N \\ t_N \end{bmatrix} + \begin{bmatrix} f_1 \\ f_2 \\ \vdots \\ f_{2N-1} \\ f_{2N} \end{bmatrix} \quad [21] $$

The new set of master coordinates is defined as:

$$ \begin{cases} s_k = a \cos \phi \\ t_k = a \sin \phi \end{cases} \quad [22] $$

where $s_k$ and $t_k$ are one pair of real linear coordinates defined in Equation (18), chosen as one of the oscillating modes.

The equations of coherence are:

with nonlinear forces:

$$ \{f\} = [Z']^{-1} \begin{Bmatrix} \{O\} \\ \{F_{NL}(\{x(\eta, \dot{\eta})\}, \{\dot{x}(\eta, \dot{\eta})\})\} \end{Bmatrix} $$

$$ \begin{cases} \dot{s}_k = s_k \Leftrightarrow \dot{a} \cos \phi - a\dot{\phi} \sin \phi \\ = \alpha_{2k-1}a \cos \phi + \beta_{2k-1}a \sin \phi + f_{2k-1} \\ \dot{t}_k = t_k \Leftrightarrow \dot{a} \sin \phi + a\dot{\phi} \cos \phi \\ = \beta_{2k}a \cos \phi + \alpha_{2k}a \sin \phi + f_{2k} \end{cases} \quad [23] $$
<figure><img src="image_12.png" /><figcaption>FIGURE 9
Fourth linear mode shape at $\Omega = 100$ rad/s: $\omega_4 = 775.684$, $\zeta_4 = 0.0561$.</figcaption></figure>

Equation (23) yields:

$$ egin{cases} 
\dot{a} = \alpha_{2k-1} a \cos^2 \phi + \alpha_{2k} a \sin^2 \phi + f_{2k-1} \cos \phi + f_{2k} \sin \phi \\ 
\dot{\phi} = -\beta_{2k-1} + (\alpha_{2k} - \alpha_{2k-1}) \sin \phi \cos \phi - \frac{f_{2k-1} \sin \phi}{a} \\ 
+ \frac{f_{2k} \cos \phi}{a} 
\end{cases} 
\quad [24] $$

The time-independent partial differential equations in $P_i$ and $Q_i$ governing the invariant manifold geometry can be written by combining Equations (21) and (25):

$$ egin{cases} 
\dot{s}_i = \frac{\partial P_i}{\partial a} \dot{a} + \frac{\partial P_i}{\partial \phi} \dot{\phi} = \alpha_{2i-1} s_i + \beta_{2i-1} t_i + f_{2i-1} & i = 1, \dots, N \\ 
\dot{t}_i = \frac{\partial Q_i}{\partial a} \dot{a} + \frac{\partial Q_i}{\partial \phi} \dot{\phi} = \beta_{2i} s_i + \alpha_{2i} t_i + f_{2i} & i \neq k 
\end{cases} 
\quad [26] $$

Following the definition of the invariant manifold, all the other d.o.f. are constrained in the $(a, \phi)$ coordinates:

$$ egin{cases} 
s_i = P_i(a, \phi) \\ 
t_i = Q_i(a, \phi) 
\end{cases} 
\quad i = 1, \dots, N \quad i \neq k 
\quad [25] $$

In those equations, $\dot{a}$ and $\dot{\phi}$ have to be replaced by expressions (24). A computation is necessary to approximate the manifold. As $P_i$ and $Q_i$ must be periodic in $\phi$ due to the transformed coordinates, they are expanded as a double series in the amplitude
<figure><img src="image_13.png" /><figcaption>FIGURE 10</figcaption></figure>

Trajectory of one end of the shaft for a very small initial disturbance and Ω = 100 rad/s.

<figure><img src="image_14.png" /><figcaption>FIGURE 11</figcaption></figure>

Trajectory of one end of the shaft for a very small initial disturbance and Ω = 2550 rad/s.
<figure><img src="image_15.png" /><figcaption>FIGURE 12
Internal resonance between mode 1 and mode 2: imaginary and real parts of the eigenvalues with $\Omega$ as a parameter.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13
First nonlinear normal mode invariant manifold: $P_4$ and $Q_4$ versus amplitude $a$ and phase $\phi$.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 14</figcaption></figure>

Motion on the first nonlinear normal mode: $t_4$ versus $s_1$ and $t_1$ obtained using both the NNM reduced-order model and the full system.

<figure><img src="image_18.png" /><figcaption>FIGURE 15</figcaption></figure>

Sixth nonlinear normal mode invariant manifold: $P_4$ and $Q_4$ versus amplitude $a$ and phase $\phi$.
and the phase:

$$ egin{cases} P_i(a, \phi) = \sum_{l=1}^{N_a} \sum_{m=1}^{N_\phi} C_i^{l,m} T_{l,m}(a, \phi) \\ Q_i(a, \phi) = \sum_{l=1}^{N_a} \sum_{m=1}^{N_\phi} D_i^{l,m} U_{l,m}(a, \phi) \end{cases} \quad i = 1, \dots, N \quad i \neq k \qquad [27] $$

where $C_i^{l,m}$ and $D_i^{l,m}$ are the unknown expansion coefficients and $T_{l,m}$ and $U_{l,m}$ are the known shape functions, typically products of functions of a over $[0, a_{max}]$ and $\phi$ over $[0, 2\pi]$. These functions are harmonic cosinus and sinus functions in $\phi$ and piecewise first-degree polynomials in a. $N_a$ and $N_\phi$ are the number of expansion functions used in the a and $\phi$ directions, respectively. The more of these functions that are used, the more accurate the approximation is. Expansions (27) are substituted into Equation (26) and a Galerkin projection (multiplication by test function and integration over the chosen domain, i.e., $[0, a_{max}] \times [0, 2\pi]$) is carried out. This yields a set of nonlinear algebraic equations in the coefficients C and D which, when solved, produces a solution that minimizes the error for the given expansion in a least squares sense.

## Results and Discussion

The above single-NNM formulation cannot handle systems with internal resonances, for which a multi-nonlinear manifold formulation would be needed. Thus, it is necessary to determine under what conditions an internal resonance can occur for the system.

An internal resonance can occur if there exists positive or negative integers $m_1, m_2, \dots, m_n$ such that $m_1\omega_1 + m_2\omega_2 + \dots + m_n\omega_n = 0$. In this case, strong nonlinear interactions between the corresponding modes may take place. To check for the presence of internal resonances in the shaft-bearings system, the loci of the eigenvalues of the linearized system are plotted in the complex plane, with the angular speed $\Omega$ as a parameter, as shown in Figure 12.

Here $\Omega$ varies from 0 to 4000 rad/s, and four distinct branches of complex conjugate pairs of oscillatory eigenvalues are observed (the other modes fall outside the scale of the plot). In order to avoid a one-to-one internal resonance, the eigenvalues must not be close to one another. Note that such a resonance precisely takes place between mode 1 and 2 for $\Omega$ approximately between 0 and 200 rad/s. It was found that in this region, the single-mode manifolds could not be calculated, and a two-mode manifold would be needed to capture the internal resonance. Hence this region was avoided.

<figure><img src="image_19.png" /><figcaption>FIGURE 16</figcaption></figure>

Motion on the sixth nonlinear normal mode: $t_4$ versus time obtained using both the NNM reduced-order model and the full system.
In order to avoid the above internal resonance, nonlinear manifolds were calculated for $\Omega = 300$ rad/s. Figure 13 shows components of the invariant manifold for the first nonlinear oscillating normal mode. Namely, the contributions $P_4$ and $Q_4$ of the fourth linear mode are depicted in terms of the amplitude and phase $\phi$. For a motion initiated on this invariant manifold, Figure 14 shows the nearly perfect agreement between the time histories calculated from the NNM single-d.o.f. reduced-order model and the full shaft bearings modal, Equation (5). Note that the latter is another possible view of a manifold, directly in the $(s,t)$-space instead of the $(a,\phi)$-space. It is important to notice that using the new NNM formulation and the real linear transformation, nonlinear manifolds can be calculated for both overdamped and oscillatory normal modes. In every computation, the nonlinear manifolds have been constructed with $N_a = N_\phi = 20$, and the dynamics on these manifolds are compared to the exact solution, using a fourth order Runge-Kutta method time-integration of these first-order differential equations with the same initial conditions.

Figures 15a and 15b depict the contributions $P_4$ and $Q_4$ of the fourth linear mode to the sixth nonlinear overdamped mode. Differences between the NNM reduced-order model and the exact solution are negligible as illustrated in Figure 16, for a motion initiated on the sixth NNM manifold.

## CONCLUSIONS

The nonlinear normal mode invariant manifold approach has been successfully extended to describe the vibration behavior of nonlinear mechanical systems with general damping, gyroscopic, and stiffness matrices. This new formulation is based on a real linear transformation and a new set of coordinates in amplitude and phase. The nonlinear Galerkin method used to solve for the invariant manifold geometry gives very good results and is quite fast. The methodology has been applied successfully to the construction of the NNM of a rotating shaft with nonlinear bearings.

Although similar results could be obtained using classical methods for nonlinear systems (e.g., harmonic balance, multiple scales), the present approach produces systematically very accurate and minimal NNM reduced-order models. It holds significant promise for the study of realistic shaft models with mass imbalance, physical phenomena such as oil whirl, and internal resonances using multi-mode manifolds.

## NOMENCLATURE

<table><tr><td>A</td><td>area of the cross section of the shaft, $2.438 \cdot 10^{-3} \text{ m}^2$</td></tr><tr><td>D1</td><td>inner diameter of the shaft, $2 \cdot 10^{-2} \text{ m}$</td></tr><tr><td>D2</td><td>outer diameter of the shaft, $5.92 \cdot 10^{-2} \text{ m}$</td></tr><tr><td>D3</td><td>inner diameter of the bearing, $5.92 \cdot 10^{-2} \text{ m}$</td></tr><tr><td>c</td><td>nominal clearance, $5 \cdot 10^{-6} \text{ m}$</td></tr><tr><td>e</td><td>eccentricity</td></tr><tr><td>E</td><td>young modulus, $2.1 \cdot 10^{11} \text{ Pa}$</td></tr></table>

<table><tr><td>ε = e/c</td><td>dimensionless eccentricity</td></tr><tr><td>{η}</td><td>vector defined by {s<sub>1</sub>,t<sub>1</sub>,...,s<sub>N</sub>,t<sub>N</sub>}</td></tr><tr><td>f</td><td>load (weight of the shaft), 186.42 N</td></tr><tr><td>{f}</td><td>vector of the bearing nonlinear forces in the (s,t)-space</td></tr><tr><td>{F<sub>NL</sub>}</td><td>vector of the bearing nonlinear forces</td></tr><tr><td>h</td><td>oil film thickness</td></tr><tr><td>I</td><td>second moment of inertia of the cross section, $1.19 \cdot 10^{-6} \text{ m}^4$</td></tr><tr><td>L</td><td>length of the shaft, 1 m</td></tr><tr><td>L<sub>b</sub></td><td>length of the bearing, $2.86 \cdot 10^{-2} \text{ m}$</td></tr><tr><td>N = Ω/2π</td><td>rotational frequency of the shaft</td></tr><tr><td>Ω</td><td>angular velocity of the shaft, rad/s</td></tr><tr><td>p</td><td>pressure in the oil film</td></tr><tr><td>φ</td><td>attitude angle</td></tr><tr><td>R = D2/2</td><td>outer radius of the shaft, $2.96 \cdot 10^{-2} \text{ m}$</td></tr><tr><td>ρ</td><td>mass density of the shaft, $7800 \text{ kg} \cdot \text{m}^{-3}$</td></tr><tr><td>μ</td><td>viscosity, $6.8 \cdot 10^{-3} \text{ N} \cdot \text{s}/\text{m}^3$</td></tr><tr><td>u, w</td><td>displacement of the shaft in X and Z direction, respectively</td></tr><tr><td>x<sub>j</sub>, z<sub>j</sub></td><td>coordinates of the center of the shaft in the bearing</td></tr><tr><td>{x}</td><td>generalized coordinates vector</td></tr><tr><td>[M<sub>1</sub>], [G<sub>1</sub>] and [K<sub>1</sub>]</td><td>mass, gyroscopic, and stiffness matrices of the shaft, respectively</td></tr><tr><td>[K<sub>2</sub>], [D<sub>2</sub>]</td><td>stiffness and damping matrices obtained by linearizing the nonlinear forces created by the oil film, respectively</td></tr></table>

## REFERENCES

* Boivin, N. 1995. Nonlinear modal analysis of structural systems using invariant manifolds. Ph.D. dissertation, Department of Mechanical Engineering, The University of Michigan.

* Craig, R. R., and Bampton, M. C. 1968. Coupling of substructures for dynamic analyses. *AIAA Journal* 6(7):1313–1319.

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

* Hirsch, M. W., and Smale, S. 1974. *Differential Equations, Dynamical Systems, and Linear Algebra*. Academic Press, New York.

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

* Kim, Y. B., and Noah, S. T. 1991b². Response and bifurcation analysis of MDOF rotor system with a strong local nonlinearity. *Nonlinear Dynamics* 2:215–234.

* Lee, C. W. 1993. *Vibration Analysis of Rotors*. Kluwer Academic Publishers, Dordrecht, Boston.

* Muszynska, A. 1986. Whirl and whip rotor/bearing stability problems, *Journal of Sound and Vibration* 110:443–462.

* Nayfeh, A. H., and Nayfeh, S. A. 1994. On nonlinear modes continuous systems. *ASME Journal of Vibration and Acoustics* 116:129–136.

* Nayfeh, A. H., and Mook, D. T. 1979. *Nonlinear Oscillations*. John Wiley & Sons, New York.
Ocvirck, F. 1952. Short bearing approximation for full journal bearings.
*NACA TN 20808*.

Pesheck, E. 2000. Reduced order modeling of nonlinear structural sys-
tems using nonlinear normal modes and invariant manifolds, Ph.D.
dissertation, Department of Mechanical Engineering, The University
of Michigan.

Rosenberg, R. M. 1966. On nonlinear vibration of systems with many
degrees of freedom. *Advances in Applied Mechanics* 9:155-242.

Shaw, S. W., and Pierre, C. 1993. Nonlinear normal modes and invariant manifolds. *Journal of Sound and Vibration* 150:170-173.

Vance, J. M. 1988. *Rotordynamics of Turbomachinery*. John Wiley &
Sons, New York.

Yamauchi, S. 1983. The nonlinear vibration of flexible rotors, first re-
port, development of a new analysis technique. *JSME 49(446) Series
C, 1862-1868.*