Lumped parameters models for the stability analysis of rotors supported on gas bearings

Abstract

In this paper different lumped parameters models are presented for the stability analysis of rotors supported on gas bearings. The analysis is carried out taking into account the damping and stiffness coefficients of journal bearings, calculated with the perturbation method. Lumped parameters models of different complexity are discussed for the description of the main rigid modes of the spindle. In case of symmetric systems, it is shown that the complexity of the model can be reduced by halving the number of the degrees of freedom. The analogy with the single mass approach is demonstrated. The considerations discussed in the paper are preliminary for the stability analysis of more complex systems, such as the ones with non-fixed bushes.

1 Introduction

Gas bearings systems suffer from the so-called half speed whirl, which is a self-excited phenomenon which occurs at high rotational speeds [1]. It involves the increase of the radius of the rotor orbits until a possible crash of the rotor with the bearings. The onset of the half-speed whirl occurs when the rotational speed exceeds the stability threshold. This phenomenon can be studied with two alternative approaches: the so-called orbit method and the linearized coefficients method.

The first method solves the time dependent Reynolds equation (RE) coupled with the equations of motion of the rotor [2]. It takes into account the non-linear effects in bearings, but it is time expensive and it is unsuitable for creating stability maps. Explicit or implicit schemes can be employed to solve the time marching problem. A noteworthy method is the Alternating Direction Implicit (ADI) scheme [3], which involves the solution of a tridiagonal system and is less costly than the implicit scheme which involves a penta-diagonal system.

The second method is based on the linearization of the bearings forces around a given operational point [4] and on the solution of the eigenproblem resulting from the rotor equations of motion [5, 6]. The linearized stiffness and damping coefficients are computed with the so-called perturbation method [7]. The improved version of the method, presented by the same author in [8], is useful for rotor-dynamics calculation of unbalance response and stability. The importance of this method was confirmed by other authors, which implemented it to carry out the stability analysis of gas bearings systems. In [9] it was employed to analyze the stability of a foil bearing and obtain stability maps, which typically show the stability threshold curve on a plot where a mass parameter is represented vs the bearing number. The classical perturbation method and an extended perturbation procedure incorporating the foil compliance explicitly were developed in [10]. The linearization method was used in [11] to compute the gas film coefficients for a high speed cryogenic turboexpander and in [12] for a heat pump turbocompressor. An interesting analysis of tilting pads is given in [13]. Many other works on the topic can be found in literature, such as [14,15,16,17]. The linearized coefficients method is more convenient from the computational point of view than the time advancing methods as the time needed to compute the linearized coefficients and evaluate the stability is much smaller than the time needed to solve the time dependent problem. Moreover, it is capable of providing a better knowledge about the underlying sources of instability, if compared with the orbit method. Stability charts can be easily plotted to show quantitatively how much the system approaches the stability threshold [6, 7, 18,19,20].

The stability evaluation is in most of cases carried out considering the so-called single-mass approach [21,22,23], which simplifies the complexity of the full rotor-bearings model and gives a first approximation of the stability threshold. Anyway, it is not evident in which cases this approach is more convenient and in which cases the full rotor model is to be preferred. To clarify this aspect, this paper shows the equivalence of the single mass approach with the simplified rotor models with 2 degrees of freedom (DOFS) and put in evidence in which conditions these models are no more accurate in predicting the stability threshold of the rotor-bearings system. In particular, it demonstrates that in case the bearings are not symmetrical with respect to the rotor center of mass, the single mass approach gives just a first estimation of the instability threshold. The error in the estimation of the onset speed increases the more unsymmetrical is the system.

In the introduction the linearized coefficients method to evaluate the stability of rotors supported on gas bearings is discussed in this paper. The full rotor model with four degrees-of-freedom (4 DOFs) system is presented, showing in which cases it is possible to reduce the number of DOFs and consider separately the cylindrical and the conical whirl motions. The 2 DOFs models are also compared with the single point-mass approach and the analogy between the two is demonstrated. Such considerations are preliminary for the stability analysis of more complex systems, such as the ones with non-fixed bushes.

1.1 The lumped parameters models

The dynamic equations of motion of the spindle considered rigid are considered in case of fixed bushings. The bearings forces are linearized around the central shaft positions and the coefficients of stiffness and damping are collected in matrixes [K] and [C]. The system dynamics is described by Eq. (1)

$$\left[M\right]\ddot{q}+\left(\left[G\right]+\left[C\right]\right)\dot{q}+\left[K\right]q={F}_{unb}$$

(1)

where [M] and [G] are the stiffness and the gyroscopic matrixes respectively, q is the state vector of the generalized displacements and F_unb is the vector of the actions due to the residual unbalance of the shaft.

Matrixes [K] and [C] depend on the linearized coefficients of the radial air films in journal bearings, which are frequency dependent.

1.2 Rigid rotor with 4 DOFs

The rigid rotor model with 4 DOFs involves all the radial degrees of freedom of the shaft. The axial degree of freedom is not considered as the axial mode and the radial modes can be considered decoupled in a first approximation; the thrust bearing is supposed to have a little influence on the whirl stability of the electro-spindle due to its small radius compared with the axial distance of the journal bearings. It is possible to consider two alternative versions of the generalized displacements vector: vector q involves both the radial displacement of the center of mass of the shaft and the rotations of the shaft axis around axes x and y:

$$q=\left[\begin{array}{cc}\begin{array}{cc}{x}_{C}& {y}_{C}\end{array}& \begin{array}{cc}{\vartheta }_{y}& -{\vartheta }_{x}\end{array}\end{array}\right]^{\prime}$$

Vector q′ involves the radial displacements of the shaft measured in the middle radial plane of the journal bearings:

$$q^{\prime}=\left[\begin{array}{cc}\begin{array}{cc}{x}_{C1}& {y}_{C1}\end{array}& \begin{array}{cc}{x}_{C2}& {y}_{C2}\end{array}\end{array}\right]^{\prime}$$

The following kinematic relationship is valid

$${q}^{\prime}=\left[\begin{array}{cc}\begin{array}{cc}1&\quad 0\\ 0&\quad 1\end{array}& \begin{array}{cc}{l}_{1}&\quad 0\\ 0&\quad {l}_{1}\end{array}\\ \begin{array}{cc}1&\quad 0\\ 0&\quad 1\end{array}& \begin{array}{cc}-{l}_{2}&\quad 0\\ 0&\quad -{l}_{2}\end{array}\end{array}\right]q$$

(2)

where l₁ and l₂ are the axial distances between the shaft center of mass and the middle plane of journal bearings, see Fig. 1. Notice that in general the system is not symmetric and distances l₁ and l₂ do not coincide.

Fig. 1

Sketch of the system’s geometry

In case q is used, the system is described by Eq. (1), while in case q’ is used, by Eq. (3)

$$\left[M^{\prime}\right]\ddot{q^{\prime}}+\left(\left[G^{\prime}\right]+\left[C^{\prime}\right]\right)\dot{q^{\prime}}+\left[K^{\prime}\right]q^{\prime}={F}_{unb}$$

(3)

The expressions of the matrixes and of the unbalance vector are reported in Appendix A1. It is worth noticing that in Eq. (1) the mass matrix is diagonal, while stiffness and damping matrixes are full. This means that in general the four equations are coupled. Similarly for system (2), in which also the mass matrix is not diagonal.

1.3 Rigid rotor with 2 DOFs, cylindrical mode

This model involves just the radial translation of the shaft, neglecting the tilting motions around axes x and y. The state vector is

$$q=\left[\begin{array}{cc}{x}_{C}& {y}_{C}\end{array}\right]^{\prime}$$

The equations of motion in matrix form are

$$\left[\begin{array}{cc}m& 0\\ 0& m\end{array}\right]\ddot{q}+\left[\begin{array}{cc}{c}_{xx1}+{c}_{xx2}& {c}_{xy1}+{c}_{xy2}\\ {c}_{yx1}+{c}_{yx2}& {c}_{yy1}+{c}_{yy2}\end{array}\right]\dot{q}+\left[\begin{array}{cc}{k}_{xx1}+{k}_{xx2}& {k}_{xy1}+{k}_{xy2}\\ {k}_{yx1}+{k}_{yx2}& {k}_{yy1}+{k}_{yy2}\end{array}\right]q=me{\omega }^{2}\left[\begin{array}{c}{\text{cos}}\,\omega t\\ {\text{sin}}\,\omega t\end{array}\right]$$

(4)

where e is the radial distance between the shaft center of mass and the center of journals (residual static unbalance), m is the mass of the shaft and k_ij and c_ij the linearized coefficients of the air films.

1.4 Rigid rotor with 2 DOFs, conical mode

In this case only the rotation of the shaft axis around x and y fixed axes are taken into account, neglecting the radial motion of the center of mass of the shaft. The state vector is

$$q=\left[\begin{array}{cc}{\vartheta }_{y}& -{\vartheta }_{x}\end{array}\right]$$

The equations of motion in matrix form are

$$\begin{aligned} \left[ {\begin{array}{*{20}c} I & 0 \\ 0 & I \\ \end{array} } \right]\ddot{q} & + \left( {\left[ {\begin{array}{*{20}c} 0 & {I_{p} \omega } \\ { - I_{p} \omega } & 0 \\ \end{array} } \right] + \left[ {\begin{array}{*{20}c} {c_{{xx1}} l_{1}^{2} + c_{{xx2}} l_{2}^{2} } & {c_{{xy1}} l_{1}^{2} + c_{{xy2}} l_{2}^{2} } \\ {c_{{yx1}} l_{1}^{2} + c_{{yx2}} l_{2}^{2} } & {c_{{yy1}} l_{1}^{2} + c_{{yy2}} l_{2}^{2} } \\ \end{array} } \right]} \right)\dot{q} \\ & + \left[ {\begin{array}{*{20}c} {k_{{xx1}} l_{1}^{2} + k_{{xx2}} l_{2}^{2} } & {k_{{xy1}} l_{1}^{2} + k_{{xy2}} l_{2}^{2} } \\ {k_{{yx1}} l_{1}^{2} + k_{{yx2}} l_{2}^{2} } & {k_{{yy1}} l_{1}^{2} + k_{{yy2}} l_{2}^{2} } \\ \end{array} } \right]q = \omega ^{2} \left[ {\begin{array}{*{20}c} {\left( {I - I_{p} } \right)\gamma \cos \omega t} \\ {\left( {I - I_{p} } \right)\gamma \sin \omega t} \\ \end{array} } \right] \\ \end{aligned}$$

(5)

where I and I_p are the transversal and the polar moments of inertia of the shaft and γ is the angle between the principal axis of inertia and the rotational axis of the shaft (residual dynamic unbalance).

1.5 Computation of gas film coefficients

The stiffness and damping coefficients of the air films are computed with the so-called perturbation method. Starting from the static pressure distribution p₀, obtained solving the discretized Reynolds equation with finite difference technique, the dynamic pressure distribution Δp is obtained as a result of the sinusoidal perturbation Δx imposed to the shaft. As sketched in Fig. 2, the coefficients are both function of the rotational speed ω and the perturbation frequency ν.

Fig. 2

Layout of the input and output variables involved in the static and in the perturbed problems

The computation of such coefficients involves the solution of a linear system, resulting after discretization of the Reynolds equation in perturbed form (6):

$$\begin{aligned} \frac{1}{{12\mu RT}}\left[ {\oint\limits_{\Gamma } {p_{0} h_{0}^{3} \vec{\nabla }\Delta p \cdot \vec{n}~d\Gamma + } \oint\limits_{\Gamma } {h_{0}^{3} \vec{\nabla }p_{0} \cdot \vec{n}\Delta pd\Gamma - } \oint\limits_{\Gamma } {3h_{0}^{2} p_{0} \Delta x\cos \vartheta \vec{\nabla }p_{0} \cdot \vec{n}~d\Gamma } } \right] \\ & \quad - \frac{1}{{RT}}\oint\limits_{\Gamma } {\left( { - p_{0} \Delta x\cos \vartheta + h_{0} \Delta p} \right)} \left[ {\begin{array}{*{20}c} {\frac{{\omega r}}{2}} \\ 0 \\ \end{array} } \right] \cdot \vec{n}~d\Gamma + G_{p}^{'} \Delta p + G_{x}^{'} \Delta x \\ & \quad = \frac{{j\nu }}{{RT}}\iint\limits_{\Sigma } {\left( {h_{0} \Delta p - p_{0} \Delta x\cos \vartheta } \right)d\Sigma } \\ \end{aligned}$$

(6)

Equation (6) represents a balance of flow in the perturbed condition; also the input flow through the supply orifices is considered in the balance, taking into account its partial derivatives with respect to pressure p (${G}_{p}^{\mathrm{^{\prime}}})$ and x (${G}_{x}^{\mathrm{^{\prime}}})$.

1.6 The Eigenvalue problem and the iterative procedure

The dynamic problem is written in the state space form, resulting in a linear system of size 2n, where n is the number of DOFs:

$$\left\{\dot{x}\right\}=\left[A\right]\left\{x\right\}$$

$$\left\{\begin{array}{c}\dot{q}\\ \ddot{q}\end{array}\right\}=\left[\begin{array}{cc}0& I\\ -{M}^{-1}K& -{M}^{-1}C\end{array}\right]\left\{\begin{array}{c}q\\ \dot{q}\end{array}\right\}$$

(7)

The stability problem is faced calculating the eigenvalues of the state space matrix A. Complex eigenvalues appear in pairs, with conjugate imaginary parts:

$$\lambda ={\lambda }_{r}\pm j{\lambda }_{i}=\eta \pm j{\omega }_{d}$$

(8)

The real part represents damping, while the imaginary part the damped frequency. The natural frequencies and the damping factors of the modes are computed with:

$${\omega }_{n}=\sqrt{{\lambda }_{r}^{2}+{\lambda }_{i}^{2}}$$

$$\zeta =\frac{-{\lambda }_{r}}{\sqrt{{\lambda }_{r}^{2}+{\lambda }_{i}^{2}}}$$

(9)

As the coefficients of the air films are frequency dependent, an iterative procedure is necessary to find the convergence between the damped frequency ${\omega }_{d}$ and the perturbation frequency ν. The following algorithm is thus implemented for each mode:

• Step1: k_ij and c_ij are evaluated at a first attempt frequency ν₀

• Step 2: the eigenvalues of system (3) are calculated, obtaining the damped frequency ω_d

• Step 3: k_ij and c_ij are evaluated at the new frequency ν  = ω_d

• Step 4: in case the new frequency coincides with the damped frequency the algorithm has reached convergence, otherwise come back to step 2.

An example of the solution of the eigenvalue problem at different rotational speeds using this algorithm is given in Fig. 3; the damped frequency and the damping factors of the four modes are plotted vs the rotational speed considering the 4 DOFs model with fixed bushings. The negative damping factor of mode 4 indicates at which rotational speed the whirl becomes unstable (about ω_th = 130 krpm). It is sufficient that one mode becomes unstable to have the unstable whirl.

Fig. 3

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed; 4 DOFs model, case a)

2 Discussion

In this paragraph, it is shown how it is possible to simplify the 4 DOFs model if suitable conditions are verified. Three cases are considered:

• (a) Different journal bearings (JBs) with unsymmetrical configuration

• (b) Different JBs with symmetrical configuration

• (c) Equal JBs with symmetrical configuration

It is shown that in the third case, in case the Coriolis terms are negligible, the equations of motion of translation and of rotation are uncoupled, so that it is possible to consider separately the cylindrical and the conical whirl motions in spite of considering the complete 4 dofs model. The bearings are supplied through two series of 10 supply holes of dia. 0.15 mm located at one fourth of the axial length. The radial film thickness is 20 µm. The rotor has mass m = 230 g, transversal moment of inertia I = 453∙10^–6 kg m² and polar moment of inertia I_p = 14.9∙10^–6 kg m².

2.1 Case a)

In this case, different JBs with unsymmetrical configuration are considered, see Table 1. The 4 DOFs model results are depicted in Fig. 3. The 2 DOFs cylindrical whirl model results are shown in Fig. 4 in the left, while the 2 DOFs conical whirl model results are shown in the right. The stability threshold obtained from the models are compared in Table 2.

Table 1 geometry considered in case a)

Fig. 4

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed for case a); 2 DOFs model, cylindrical whirl in the left; 2 DOFs model, conical whirl in the right

Table 2 comparison of the stability threshold according to the considered models; case a)

By a comparison of the 2 DOFs models with the 4 DOFs model, the damped frequencies and the damping ratio of the modes do not coincide. The reason is that in this case the cylindrical and the conical whirl motions are coupled and the 2 DOFs models, which take into account only the pure cylindrical or the conical whirl, result to be inaccurate. The complete 4 DOFs model is to be preferred as it treats the coupled cylindrical and conical whirl.

2.2 Case b)

In this case, different JBs with symmetrical configuration are considered, see Table 3.

Table 3 geometry considered in case b)

The 4 DOFs model results are depicted in Fig. 5. The 2 DOFs cylindrical whirl model results are shown in Fig. 6 (left), while the 2 DOFs conical whirl model results are shown in the right. The stability threshold obtained from the models are compared in Table 4.

Fig. 5

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed; 4 DOFs model, case b)

Fig. 6

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed for case b); 2 DOFs model, cylindrical whirl in the left; 2 DOFs model, conical whirl in the right

Table 4 comparison of the stability threshold according to the considered models; case b)

In this case the curves representing the damped frequency of the modes and the curves of the damping factor are more similar if the 4 DOFs model is compared with the 2 DOFs models. Due to the increased symmetry of the system, the 2 DOFs model give a better approximation of the 4 DOFs model with respect to case a). Anyway, the results do not yet coincide.

2.3 Case c)

In this case, equal JBs with symmetrical configuration are considered, see Table 5.

Table 5 geometry considered in case c)

The 4 DOFs model results are depicted in Fig. 7. The 2 DOFs cylindrical whirl model results are shown in Fig. 8 in the left, while the 2 DOFs conical whirl model results are shown in the right.

Fig. 7

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed; 4 DOFs model, case c)

Fig. 8

Damped frequencies ω_d and damping factors $\zeta$ vs rotational speed for case c); 2 DOFs model, cylindrical whirl in the left; 2 DOFs model, conical whirl in the right

The stability threshold obtained from the models are compared in Table 6, while Table 7 compares the stiffness and damping matrixes in the three cases. Stiffness coefficients are expressed in N/m, while damping coefficients in Ns/m. The table shows how matrixes [K] and [C] become band diagonal the more symmetrical the geometry becomes.

Table 6 comparison of the stability threshold according to the considered models, case c)

Table 7 comparison of matrixes [K] and [C] in cases a), b) and c)

In this case, due to the complete symmetry of the journal bearings with respect to the center of mass of the shaft, the results of the 2 DOFs models perfectly match with the results of the 4 DOFs model. The curves representing the damped frequencies and the damping factors coincide.

The cylindrical and the conical motions are uncoupled in case the Coriolis terms are neglected and the system’s geometry satisfies condition (10) or condition (11). In the first case there is perfect symmetry: distances l₁ and l₂ must coincide and the journal bearings must have the same geometry and the same radial film thickness (same air film coefficients):

$$\left\{\begin{array}{l}\begin{array}{l}{l}_{1}={l}_{2}\\ {k}_{ij, 1}={k}_{ij,2} \end{array}\\ {c}_{ij, 1}={c}_{ij,2}\end{array}\right.$$

(10)

Condition (11) is more general, as the center of mass of the shaft could be not in the middle of bearings.

$$\left\{\begin{array}{c}{k}_{ij, 1}{l}_{1}={k}_{ij,2}{l}_{2}\\ {c}_{ij, 1}{l}_{1}={c}_{ij,2}{l}_{2}\end{array}\right.$$

(11)

This is evident from the stiffness and damping matrixes of the 4dofs model (see appendix A1) which reduce to band diagonal matrixes. For instance, in condition (10) [K] is reduced into this form:

$$\left[ {\begin{array}{*{20}c} {2k_{{xx}} } & {2k_{{xy}} } & 0 & 0 \\ {2k_{{yx}} } & {2k_{{yy}} } & 0 & 0 \\ 0 & 0 & {2k_{{xx}} l^{2} } & {2k_{{xy}} l^{2} } \\ 0 & 0 & {2k_{{yx}} l^{2} } & {2k_{{yy}} l^{2} } \\ \end{array} } \right]$$

(12)

It is thus demonstrated in which cases the full rotor model can be reduced: in case of complete symmetry (condition 10) or more in general according to condition (11). Of course, in case the system does not exactly satisfy such conditions but the asymmetry is small or condition (11) is almost satisfied, the two modes can be still studied separately, being the out-of-diagonal elements negligible. In the other cases, the full model is to be preferred as it better estimates the stability threshold.

2.4 The single mass approach

In literature the so-called single mass approach is known, which allows a quick estimation of the stability of rotor-bearings systems. As explained in [23, 24], the approach is quite common to perform the stability analysis of a rotor-journal bearings systems. According to this approach, the shaft mass is distributed in correspondence of the journal bearings considering a cylindrical motion or a conical motion. The masses m_s1 and m_s2 change in case the cylindrical or the conical whirl are considered, see Fig. 9.

Fig. 9

Lumped masses in correspondence of journal bearings

Such expressions are compared with the masses that appear in [M’] in case a pure cylindrical or a pure conical motion are considered, see Appendix A2 for the detailed passages. It is thus demonstrated the analogy of the single mass approach with the simplified 2 DOFs models, in case of proper symmetry of the system.

3 Conclusions

This paper illustrates in details the linearized stability analysis of rotor-gas bearings systems. It shows in which cases a 4dofs rotor-bearings model can be simplified in a 2dofs model considering alternatively the cylindrical whirl or the conical whirl. In these cases the stability analysis can be carried out on each separate bearing concentrating on it a fraction of the shaft mass. The analogy of the reduced models with the so-called single mass approach is demonstrated. The models can be further extended to the stability analysis of gas bearings systems with non-fixed bushings.

Appendices

Appendix A1: 4 DOFS MODEL

The matrixes of the 4 DOFs model of the rigid shaft are here explicitated. In case the generalized displacement vector q is used, the following expressions result:

$$\begin{array}{*{20}l} {\left[ M \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} m & 0 \\ 0 & m \\ \end{array} } & {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} I & 0 \\ 0 & I \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]} \hfill & {\left[ K \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {k_{{xx1}} + k_{{xx2}} } & {k_{{xy1}} + k_{{xy2}} } \\ {k_{{yx1}} + k_{{yx2}} } & {k_{{yy1}} + k_{{yy2}} } \\ \end{array} } & {\begin{array}{*{20}c} {k_{{xx1}} l_{1} - k_{{xx2}} l_{2} } & {k_{{xy1}} l_{1} - k_{{xy2}} l_{2} } \\ {k_{{yx1}} l_{1} - k_{{yx2}} l_{2} } & {k_{{yy1}} l_{1} - k_{{yy2}} l_{2} } \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {k_{{xx1}} l_{1} - k_{{xx2}} l_{2} } & {k_{{xy1}} l_{1} - k_{{xy2}} l_{2} } \\ {k_{{yx1}} l_{1} - k_{{yx2}} l_{2} } & {k_{{yy1}} l_{1} - k_{{yy2}} l_{2} } \\ \end{array} } & {\begin{array}{*{20}c} {k_{{xx1}} l_{1}^{2} + k_{{xx2}} l_{2}^{2} } & {k_{{xy1}} l_{1}^{2} + k_{{xy2}} l_{2}^{2} } \\ {k_{{yx1}} l_{1}^{2} + k_{{yx2}} l_{2}^{2} } & {k_{{yy1}} l_{1}^{2} + k_{{yy2}} l_{2}^{2} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]} \hfill \\ {\left[ G \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & {I_{p} \omega } \\ { - I_{p} \omega } & 0 \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]} \hfill & {\left[ C \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} {c_{{xx1}} + c_{{xx2}} } & {c_{{xy1}} + c_{{xy2}} } \\ {c_{{yx1}} + c_{{yx2}} } & {c_{{yy1}} + c_{{yy2}} } \\ \end{array} } & {\begin{array}{*{20}c} {c_{{xx1}} l_{1} - c_{{xx2}} l_{2} } & {c_{{xy1}} l_{1} - c_{{xy2}} l_{2} } \\ {c_{{yx1}} l_{1} - c_{{yx2}} l_{2} } & {c_{{yy1}} l_{1} - c_{{yy2}} l_{2} } \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {c_{{xx1}} l_{1} - c_{{xx2}} l_{2} } & {c_{{xy1}} l_{1} - c_{{xy2}} l_{2} } \\ {c_{{yx1}} l_{1} - c_{{yx2}} l_{2} } & {c_{{yy1}} l_{1} - c_{{yy2}} l_{2} } \\ \end{array} } & {\begin{array}{*{20}c} {c_{{xx1}} l_{1}^{2} + c_{{xx2}} l_{2}^{2} } & {c_{{xy1}} l_{1}^{2} + c_{{xy2}} l_{2}^{2} } \\ {c_{{yx1}} l_{1}^{2} + c_{{yx2}} l_{2}^{2} } & {c_{{yy1}} l_{1}^{2} + c_{{yy2}} l_{2}^{2} } \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]} \hfill \\ \end{array}$$

(13)

$${F}_{unb}={\omega }^{2}\left[\begin{array}{c}\begin{array}{c}me{\text{cos}}\,\omega t\\ me{\text{sin}}\,\omega t\end{array}\\ \begin{array}{c}\left(I-{I}_{p}\right)\gamma {\text{cos}}\left(\omega t+{\varphi }_{1}\right)\\ \left(I-{I}_{p}\right)\gamma {\text{sin}}\left(\omega t+{\varphi }_{1}\right)\end{array}\end{array}\right]$$

In case the generalized displacement vector q’ is used, the following expressions result:

$$\begin{array}{*{20}l} {\left[ {M^\prime } \right] = \frac{1}{{l_{1} + l_{2} }}\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {ml_{2} } & 0 \\ 0 & {ml_{2} } \\ \end{array} } & {\begin{array}{*{20}c} {ml_{1} } & 0 \\ 0 & {ml_{1} } \\ \end{array} } \\ {\begin{array}{*{20}c} I & 0 \\ 0 & I \\ \end{array} } & {\begin{array}{*{20}c} { - I} & 0 \\ 0 & { - I} \\ \end{array} } \\ \end{array} } \right]} \hfill & {\left[ {K^\prime } \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {k_{{xx1}} } & {k_{{xy1}} } \\ {k_{{yx1}} } & {k_{{yy1}} } \\ \end{array} } & {\begin{array}{*{20}c} {k_{{xx2}} } & {k_{{xy2}} } \\ {k_{{yx2}} } & {k_{{yy2}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {k_{{xx1}} l_{1} } & {k_{{xy1}} l_{1} } \\ {k_{{yx1}} l_{1} } & {k_{{yy1}} l_{1} } \\ \end{array} } & {\begin{array}{*{20}c} { - k_{{xx2}} l_{2} } & { - k_{{xy2}} l_{2} } \\ { - k_{{yx2}} l_{2} } & { - k_{{yy2}} l_{2} } \\ \end{array} } \\ \end{array} } \right]} \hfill \\ {\left[ {G^\prime } \right] = \frac{1}{{l_{1} + l_{2} }}\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & {I_{p} \omega } \\ { - I_{p} \omega } & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & { - I_{p} \omega } \\ {I_{p} \omega } & 0 \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]} \hfill & {\left[ {C^\prime } \right] = \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {c_{{xx1}} } & {c_{{xy1}} } \\ {c_{{yx1}} } & {c_{{yy1}} } \\ \end{array} } & {\begin{array}{*{20}c} {c_{{xx2}} } & {c_{{xy2}} } \\ {c_{{yx2}} } & {c_{{yy2}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {c_{{xx1}} l_{1} } & {c_{{xy1}} l_{1} } \\ {c_{{yx1}} l_{1} } & {c_{{yy1}} l_{1} } \\ \end{array} } & {\begin{array}{*{20}c} { - c_{{xx2}} l_{2} } & { - c_{{xy2}} l_{2} } \\ { - c_{{yx2}} l_{2} } & { - c_{{yy2}} l_{2} } \\ \end{array} } \\ \end{array} } \right]} \hfill \\ \end{array}$$

(14)

Appendix A2: analogy with the single-mass approach

2.1 Cylindrical whirl

Starting from system (3) and imposing the cylindrical whirl it is.

$${x}_{C1}={x}_{C2}={x}_{C};{y}_{C1}={y}_{C2}={y}_{C};\gamma =0$$

The system becomes

$$\frac{1}{{l}_{1}+{l}_{2}}\left[\begin{array}{c}\begin{array}{cc}m{l}_{2}+m{l}_{1}& 0\\ 0& m{l}_{2}+m{l}_{1}\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_{C}\\ {\ddot{y}}_{C}\end{array}\right]+\left[\begin{array}{c}\begin{array}{cc}{k}_{xx1}+{k}_{xx2}& {k}_{xy1}+{k}_{xy2}\\ {k}_{yx1}+{k}_{yx2}& {k}_{yy1}+{k}_{yy2}\end{array}\\ \begin{array}{cc}{k}_{xx1}{l}_{1}-{k}_{xx2}{l}_{2}& {k}_{xy1}{l}_{1}-{k}_{xy2}{l}_{2}\\ {k}_{yx1}{l}_{1}-{k}_{yx2}{l}_{2}& {k}_{yy1}{l}_{1}-{k}_{yy2}{l}_{2}\end{array}\end{array}\right]\left[\begin{array}{c}{x}_{C}\\ {y}_{C}\end{array}\right]={\omega }^{2}\left[\begin{array}{c}\begin{array}{c}me{\text{cos}}\,\omega t\\ me{\text{sin}}\, \omega t\end{array}\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]$$

(15)

where for brevity the damping matrix is not indicated as it has a structure similar to the stiffness matrix. In case of symmetry (${k}_{ij1}\approx {k}_{ij2}$ and ${l}_{1}\approx {l}_{2}$, or more generally ${k}_{ij1}{l}_{1}\approx {k}_{ij2}{l}_{2}$) the last two lines are zeroed and the system becomes

$$\left[\begin{array}{cc}{m}_{s1}+{m}_{s2}& 0\\ 0& {m}_{s1}+{m}_{s2}\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_{C}\\ {\ddot{y}}_{C}\end{array}\right]+\left[\begin{array}{cc}{k}_{xx1}+{k}_{xx2}& {k}_{xy1}+{k}_{xy2}\\ {k}_{yx1}+{k}_{yx2}& {k}_{yy1}+{k}_{yy2}\end{array}\right]\left[\begin{array}{c}{x}_{C}\\ {y}_{C}\end{array}\right]={\omega }^{2}\left[\begin{array}{c}\left({m}_{s1}{e}_{1}+{m}_{s2}{e}_{2}\right){\text{cos}}\omega t\\ \left({m}_{s1}{e}_{1}+{m}_{s2}{e}_{2}\right){\text{sin}}\omega t\end{array}\right]$$

(16)

where m_s1 and m_s2 are the masses in correspondence of bearings, see Table 

Table 8 expressions of the masses concentrated in correspondence to the journal bearings

8, and e1 and e2 their unbalance. It must be

$${m}_{s1}{e}_{1}+{m}_{s2}{e}_{2}=me$$

(17)

This is the equilibrium of forces considering the two masses together. It is equivalent to the force equilibrium of the masses taken singularly, as no radial force is transmitted by the shaft; in this case, the inertia force of the single mass is supported by the correspondent bearing. The advantage of considering the mass concentrated in correspondence of bearings is that the stability of the single bearing can be evaluated separately, with a reduced number of DOFs.

$$\left[\begin{array}{cc}{m}_{si}& 0\\ 0& {m}_{si}\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_{C}\\ {\ddot{y}}_{C}\end{array}\right]+\left[\begin{array}{cc}{k}_{xxi}& {k}_{xyi}\\ {k}_{yxi}& {k}_{yyi}\end{array}\right]\left[\begin{array}{c}{x}_{C}\\ {y}_{C}\end{array}\right]={\omega }^{2}{m}_{si}{e}_{i}\left[\begin{array}{c}{\text{cos}}\omega t\\ {\text{sin}}\omega t\end{array}\right]$$

(18)

2.2 Conical whirl

Starting from system (3) and imposing the conical whirl it is.

${x}_{C1}={\vartheta }_{y}{l}_{1};{y}_{C1}=-{\vartheta }_{x}{l}_{1};{x}_{C2}=-{\vartheta }_{y}{l}_{2};{x}_{C2}=-{\vartheta }_{y}{l}_{2};{y}_{C2}={\vartheta }_{x}{l}_{2};e\hspace{0.17em}=\hspace{0.17em}0$

As $\frac{{x}_{C1}}{{l}_{1}}=-\frac{{x}_{C2}}{{l}_{2}}$ some elements in the mass matrix are zeroed:

$$\begin{aligned} \frac{1}{{l_{1} + l_{2} }}\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } \\ {\begin{array}{*{20}c} I & 0 \\ 0 & I \\ \end{array} } & {\begin{array}{*{20}c} { - I} & 0 \\ 0 & { - I} \\ \end{array} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\ddot{x}_{{C1}} } \\ {\ddot{y}_{{C1}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\ddot{x}_{{C2}} } \\ {\ddot{y}_{{C2}} } \\ \end{array} } \\ \end{array} } \right] + \frac{1}{{l_{1} + l_{2} }}\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & 0 \\ 0 & 0 \\ \end{array} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 & {I_{p} \omega } \\ { - I_{p} \omega } & 0 \\ \end{array} } & {\begin{array}{*{20}c} 0 & { - I_{p} \omega } \\ {I_{p} \omega } & 0 \\ \end{array} } \\ \end{array} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {\dot{x}_{{C1}} } \\ {\dot{y}_{{C1}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {\dot{x}_{{C2}} } \\ {\dot{y}_{{C2}} } \\ \end{array} } \\ \end{array} } \right] \\ & + \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {k_{{xx1}} } & {k_{{xy1}} } \\ {k_{{yx1}} } & {k_{{yy1}} } \\ \end{array} } & {\begin{array}{*{20}c} {k_{{xx2}} } & {k_{{xy2}} } \\ {k_{{yx2}} } & {k_{{yy2}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {k_{{xx1}} l_{1} } & {k_{{xy1}} l_{1} } \\ {k_{{yx1}} l_{1} } & {k_{{yy1}} l_{1} } \\ \end{array} } & {\begin{array}{*{20}c} { - k_{{xx2}} l_{2} } & { - k_{{xy2}} l_{2} } \\ { - k_{{yx2}} l_{2} } & { - k_{{yy2}} l_{2} } \\ \end{array} } \\ \end{array} } \right]\left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} {x_{{C1}} } \\ {y_{{C1}} } \\ \end{array} } \\ {\begin{array}{*{20}c} {x_{{C2}} } \\ {y_{{C2}} } \\ \end{array} } \\ \end{array} } \right] = \omega ^{2} \left[ {\begin{array}{*{20}c} {\begin{array}{*{20}c} 0 \\ 0 \\ \end{array} } \\ {\begin{array}{*{20}c} {\left( {I - I_{p} } \right)\gamma {\text{cos}}\left( {\omega t + \varphi _{1} } \right)} \\ {\left( {I - I_{p} } \right)\gamma {\text{sin}}\left( {\omega t + \varphi _{1} } \right)} \\ \end{array} } \\ \end{array} } \right] \\ \end{aligned}$$

(19)

In case of symmetry (${k}_{ij1}\approx {k}_{ij2}$ and ${l}_{1}\approx {l}_{2}$, or more generally ${k}_{ij1}{l}_{1}\approx {k}_{ij2}{l}_{2}$) the first two lines are zeroed and the system can be further simplified neglecting the Coriolis terms:

$$\frac{1}{{l}_{1}+{l}_{2}}\left[\begin{array}{cc}\begin{array}{cc}I& 0\\ 0& I\end{array}& \begin{array}{cc}-I& 0\\ 0& -I\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\ddot{x}}_{C1}\\ {\ddot{y}}_{C1}\end{array}\\ \begin{array}{c}{\ddot{x}}_{C2}\\ {\ddot{y}}_{C2}\end{array}\end{array}\right]+\left[\begin{array}{cc}\begin{array}{cc}{k}_{xx1}{l}_{1}& {k}_{xy1}{l}_{1}\\ {k}_{yx1}{l}_{1}& {k}_{yy1}{l}_{1}\end{array}& \begin{array}{cc}-{k}_{xx2}{l}_{2}& -{k}_{xy2}{l}_{2}\\ -{k}_{yx2}{l}_{2}& -{k}_{yy2}{l}_{2}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{x}_{C1}\\ {y}_{C1}\end{array}\\ \begin{array}{c}{x}_{C2}\\ {y}_{C2}\end{array}\end{array}\right]={\omega }^{2}\left(I-{I}_{p}\right)\gamma \left[\begin{array}{c}{\text{cos}}\left(\omega t+{\varphi }_{1}\right)\\ {\text{sin}}\left(\omega t+{\varphi }_{1}\right)\end{array}\right]$$

(20)

Substituting the expressions of masses m_s1 and m_s2 from Table 8 for conical whirl we obtain

$$\left[\begin{array}{cc}\begin{array}{cc}{m}_{s1}{l}_{1}& 0\\ 0& {m}_{s1}{l}_{1}\end{array}& \begin{array}{cc}-{m}_{s2}{l}_{2}& 0\\ 0& -{m}_{s2}{l}_{2}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{\ddot{x}}_{C1}\\ {\ddot{y}}_{C1}\end{array}\\ \begin{array}{c}{\ddot{x}}_{C2}\\ {\ddot{y}}_{C2}\end{array}\end{array}\right]+\left[\begin{array}{cc}\begin{array}{cc}{k}_{xx1}{l}_{1}& {k}_{xy1}{l}_{1}\\ {k}_{yx1}{l}_{1}& {k}_{yy1}{l}_{1}\end{array}& \begin{array}{cc}-{k}_{xx2}{l}_{2}& -{k}_{xy2}{l}_{2}\\ -{k}_{yx2}{l}_{2}& -{k}_{yy2}{l}_{2}\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{x}_{C1}\\ {y}_{C1}\end{array}\\ \begin{array}{c}{x}_{C2}\\ {y}_{C2}\end{array}\end{array}\right]=\left({m}_{s1}{e}_{1}{l}_{1}-{m}_{s2}{e}_{2}{l}_{2}\right)\left[\begin{array}{c}{\text{cos}}\left(\omega t+{\varphi }_{1}\right)\\ {\text{sin}}\left(\omega t+{\varphi }_{1}\right)\end{array}\right]$$

(21)

where it is

$$\left(I-{I}_{p}\right)\gamma =\left({m}_{s1}{e}_{1}{l}_{1}-{m}_{s2}{e}_{2}{l}_{2}\right)$$

(22)

Equations 21 represent to the equilibrium of moments of the masses concentrated in correspondence of bearings. As before, the equilibrium can be written for the single mass taken separately as no radial force is transmitted by the shaft and the inertia force of each mass is balanced by the corresponding bearing. Dividing by the arm l_i, the moment equilibrium simplifies into a force equilibrium.

$$\left[\begin{array}{cc}{m}_{si}{l}_{i}& 0\\ 0& {m}_{si}{l}_{i}\end{array}\right]\left[\begin{array}{c}{\ddot{x}}_{Ci}\\ {\ddot{y}}_{Ci}\end{array}\right]+\left[\begin{array}{cc}{k}_{xxi}{l}_{i}& {k}_{xyi}{l}_{i}\\ {k}_{yxi}{l}_{i}& {k}_{yyi}{l}_{i}\end{array}\right]\left[\begin{array}{c}{x}_{Ci}\\ {y}_{Ci}\end{array}\right]={m}_{s1}{e}_{1}{l}_{i}\left[\begin{array}{c}{\text{cos}}\left(\omega t+{\varphi }_{1}\right)\\ {\text{sin}}\left(\omega t+{\varphi }_{1}\right)\end{array}\right]$$

(23)