The effect of time-periodic base angular motions upon dynamic response of asymmetric rotor systems

Abstract

The effect of time-periodic base angular motions upon dynamic response of asymmetric rotor systems is investigated in this article. Both the disk and rotating shaft with asymmetric cross sections are considered. The rigid rotor base rotates around the three axis with time-varying angular speeds. Equations of motion for every components (asymmetric disk and shaft elements) are derived to assembe the finite element model of the whole rotor-bearing system. Due to the rotor asymmetries and angular motions of the rigid base, the finite element model has time-variable mass, gyroscopic, and stiffness coefficients. Parametric instability regions and rotor center orbits are compared with references’ results for model verifications. Then, forced responses of the system under unbalanced excitations, shaft asymmetry, and base angular motions are obtained. Effects of various types of angular motions, amplitudes, and frequency of base motions upon the response spectra and external resonances are discussed. As long as the periodic base angular motion is considered, both the response spectra and external resonances change significantly. Additional frequencies, including base frequency and combined frequencies, are found in the response spectra. Besides the critical and half-critical resonances (for the asymmetric rotor system), there are multiple resonances after considering the periodic base motions.

Introduction

It is common to find the rotor with asymmetric cross sections in the field of rotating machinery. Shaft with rectangular cross section, open-edge crack, or a key way and shaft of a two-pole turbo-generator are typical examples of asymmetric rotors. The non-symmetric cross sections cause the inertia and stiffness in two principal planes different and induce parametric excitations to the rotor system. In a remarkable paper presented in 1933, the asymmetric properties in rotor system was first noticed by Smith. ¹ In the 1960s, Crandall and Brosens ² and Yamamoto and colleagues^3,4 first conducted quantitative studies on the dynamic behaviors of asymmetric rotors. Since then, extensive work has been done on the asymmetric rotor dynamics.^5–18 According to the types of rotor models used in the analysis, current literature could be divided into three groups: the Jeffcott rotor model (or lumped parameter rotor model),^5–9 the transfer matrix model,^10–12 and the finite element (FE) rotor model (FE model (FEM)).^13–18 If one considers the research contents, the current literature could be distinguished into parametric instability analysis^5–11,15,17 and steady-state response analysis.^{12–14,16,18}

It should be noted that the above-mentioned literature all assumed that the rotor supports are fixed without motions. In reality, the rotors supported by moving bases are also encountered frequently. Rotating machines installed on transportation systems belong to such type of rotors. Over the past 20 years, the dynamic behaviors of general rotor systems under moving bases have been paid sufficient attentions by various scholars. Considering a simple Jeffcott rotor installed in an aircraft, Lin and Meng ¹⁹ first studied the large linear and angular base motions induced by the maneuver flight. They discovered that the flying status of the aircraft has significant effect on the shaft’s unbalanced response. Lee et al. ²⁰ considered a base-transferred shock force and utilized a direction integration scheme to obtain the rotor’s transient response. Driot et al. ²¹ and Duchemin et al. ²² focused their study on the whirling orbits and parametric instability caused by base angular motions. Semi-analytical expressions of instability boundary of a simple Jeffcott rotor system were obtained using the multiple scales method. Later, they also considered the random base motions for stochastic rotor dynamic analysis. ²³ El-Saeidy and Sticher ²⁴ analyzed the dynamic behaviors of a rigid rotor-bearing system under both imbalance and harmonic base excitations. Saha et al. ²⁵ calculated the dynamic response of an uniformly spinning shaft with a non-central disk mounted on a precessing base. Transverse cracks are frequently found in the rotor’s shaft. Recently, Lin et al., ²⁶ Jing et al., ²⁷ and Yang et al. ²⁸ studied the coupling effects of moving base on the cracked rotor dynamic behaviors. Through sensitivity analysis, Han and Chu ²⁹ discovered some unique spectra for detecting various transverse cracks in base-excited rotating machinery. For a rotor system with rub-impact fault and cubic stiffness supports, Hou et al. ³⁰ found some typical nonlinear phenomena induced by the aircraft hovering flight. As the base motions would induce severe vibrations to rotor system under certain conditions, several researchers also presented the active control method using the electromagnetic actuators.^31–33

From the above literature, one can see that both linear and angular base motions not only cause external force excitations but also internal excitations to the rotor-bearing system. It is well known that the rotor asymmetry results in parametric excitations to the system. Thus, it is foreseeable that the dynamic behavior of flexible rotor-bearings would be greatly changed after considering both the moving base and rotor asymmetry simultaneously. To the knowledge of the authors, no publication is available in the open literature that reports the effect of moving bases on the dynamic response of asymmetric rotor-bearing systems.

Thus, this article tries to reveal the effect of time-periodic base angular motions upon dynamic response of asymmetric rotor systems. Both the disk and rotating shaft with asymmetric cross sections are considered. The rotor base is assumed to be rigid and rotates around the three axis with time-varying angular speeds. Utilizing the energy theorem and Lagrange’s principle, the equations of motions for every components (asymmetric disk and shaft elements) are derived to assemble the FEM of the whole rotor-bearing system. Due to the rotor asymmetries and angular motions of the rigid base, the FEM has time-variable mass, gyroscopic, and stiffness coefficients. Before numerical simulations, the established rotor-bearing model is validated through instability and steady-state response analysis with reference results.^16,21 Then, forced steady-state responses of the system under unbalanced excitations, shaft asymmetry, and base angular motions are analyzed numerically. The effects of various types of angular motions, amplitudes, and frequency of base motions upon the response spectra and external resonances are discussed. Finally, some conclusions are summarized.

FEM

In order to derive the FEM, the following assumptions should be explained:

1. Rigid base is assumed and its mass center is at the left bearing (see Figure 1).

2. There are three coordinate systems defined for modeling both rotor and base motion: the inertial frame of reference $(X_0-Y_0-Z_0)$ , frame of moving base $(X_b-Y_b-Z_b)$ , and rotor frame $(X-Y-Z)$ , which are also identified in Figure 1.

3. The rotation speed of the rotor shaft is assumed to be constant and denoted by ${\Omega }_s$ . The rigid motion rotates around the three axis of the frame $X_b-Y_b-Z_b$ and the time-variable angular motions are denoted by ${\alpha }_b$ (pitching), ${\beta }_b$ (rolling), and ${\gamma }_b$ (yawing).

Figure 1.

Asymmetric rotor-bearing system on moving base with various coordinate systems.

The asymmetric disk

Four degrees of freedom, including two translations $(x_d,z_d)$ and two rotations $({\theta }_{\mathrm{dz}},{\theta }_{\mathrm{dx}})$ , are considered for the rigid asymmetric disk node. Its nodal displacement vector is then written by ${\mathbf{q}}_d=[x_d{z}_d{\theta }_{\mathrm{dz}}{\theta }_{\mathrm{dx}}{]}^T$ . It is easy to write the expressions of kinetic energy of a rigid disk as follows

$$T_d=\frac{1}{2}{m}_d({\stackrel{·}{u}}_d^2+{\stackrel{·}{v}}_d^2+{\stackrel{·}{w}}_d^2)+\frac{1}{2}(I_{\mathrm{dx}}{\omega }_{\mathrm{dx}}^2+I_p{\omega }_{\mathrm{dy}}^2+I_{\mathrm{dz}}{\omega }_{\mathrm{dz}}^2)$$

(1)

where $m_d,I_{\mathrm{dx}},I_{\mathrm{dz}},\mathrm{and}{I}_p$ are, respectively, the mass, two diametral mass-moments of inertia, and polar mass-moment of inertia of the disk located at a distance $l_d$ from the left bearing. For the asymmetric disk, one has $I_{\mathrm{dx}}\ne I_{\mathrm{dz}}$ . The mean value of the diametral moments of inertia is introduced as $I_m=1/2(I_{\mathrm{dx}}+I_{\mathrm{dz}})$ , and the relative inertia asymmetry of the disk is defined by ${\delta }_d=(I_{\mathrm{dx}}-I_{\mathrm{dz}})/2I_m$ . The inertia ratio is defined as ${\chi }_d=I_p/I_m$ . The ${\stackrel{·}{u}}_d,{\stackrel{·}{v}}_d,\mathrm{and}{\stackrel{·}{w}}_d$ and ${\omega }_{\mathrm{dx}},{\omega }_{\mathrm{dy}},\mathrm{and}{\omega }_{\mathrm{dz}}$ are, respectively, absolute translational and angular velocities of disk node with respect to the inertial frame $X_0-Y_0-Z_0$ . For the absolute translational velocities, one can write

$$\begin{array}{c}{\stackrel{·}{u}}_d={\stackrel{·}{\beta }}_b{z}_d-{\stackrel{·}{\gamma }}_b{l}_d+{\stackrel{·}{x}}_d\hfill \\ {\stackrel{·}{v}}_d={\stackrel{·}{\gamma }}_b{x}_d-{\stackrel{·}{\alpha }}_b{z}_d\hfill \\ {\stackrel{·}{w}}_d={\stackrel{·}{\alpha }}_b{l}_d-{\stackrel{·}{\beta }}_b{x}_d+{\stackrel{·}{z}}_d\hfill \end{array}$$

(2)

Three angles, ${\theta }_{\mathrm{dz}},{\theta }_{\mathrm{dx}},\mathrm{and}\Phi$ , are utilized to describe the movement of rotor frame $X-Y-Z$ with respect to the base frame $X_b-Y_b-Z_b$ . As the disk rotates constantly with angular speed ${\Omega }_s$ , one has $\Phi ={\Omega }_{s}t$ . There are three steps in order to orientate the disk. First, it is rotated by angle ${\theta }_{\mathrm{dz}}$ around axis $Z_b$ , then by angle ${\theta }_{\mathrm{dx}}$ around the new axis $X$ , and finally by angle $\Phi$ around the new axis $Y$ . Thus, the absolute angular velocities of the disk could be derived by transforming the coordinates as

$$\begin{array}{cc}\multicolumn{1}{c}{{\omega }_{\mathrm{dx}}=}& ({\stackrel{·}{\alpha }}_b\cos {\theta }_{\mathrm{dz}}+{\stackrel{·}{\beta }}_b\sin {\theta }_{\mathrm{dz}}+{\stackrel{·}{\theta }}_{\mathrm{dx}})\cos \Phi \hfill \\ \multicolumn{1}{c}{}& +(({\stackrel{·}{\alpha }}_b\sin {\theta }_{\mathrm{dz}}-{\stackrel{·}{\beta }}_b\cos {\theta }_{\mathrm{dz}})\sin {\theta }_{\mathrm{dx}}\hfill \\ \multicolumn{1}{c}{}& +(\stackrel{·}{\gamma }+{\stackrel{·}{\theta }}_{\mathrm{dz}})\cos {\theta }_{\mathrm{dx}})\sin \Phi \hfill \\ \multicolumn{1}{c}{{\omega }_{\mathrm{dy}}=}& (-{\stackrel{·}{\alpha }}_b\sin {\theta }_{\mathrm{dz}}+{\stackrel{·}{\beta }}_b\cos {\theta }_{\mathrm{dz}})\cos {\theta }_{\mathrm{dx}}\hfill \\ \multicolumn{1}{c}{}& +(\stackrel{·}{\gamma }+{\stackrel{·}{\theta }}_{\mathrm{dz}})\sin {\theta }_{\mathrm{dx}})+\stackrel{·}{\Phi }\hfill \\ \multicolumn{1}{c}{{\omega }_{\mathrm{dz}}=}& -({\stackrel{·}{\alpha }}_b\cos {\theta }_{\mathrm{dz}}+{\stackrel{·}{\beta }}_b\sin {\theta }_{\mathrm{dz}}+{\stackrel{·}{\theta }}_{\mathrm{dx}})\sin \Phi \hfill \\ \multicolumn{1}{c}{}& +(({\stackrel{·}{\alpha }}_b\sin {\theta }_{\mathrm{dz}}-{\stackrel{·}{\beta }}_b\cos {\theta }_{\mathrm{dz}})\sin {\theta }_{\mathrm{dx}}\hfill \\ \multicolumn{1}{c}{}& +(\stackrel{·}{\gamma }+{\stackrel{·}{\theta }}_{\mathrm{dz}})\cos {\theta }_{\mathrm{dx}})\cos \Phi \hfill \end{array}$$

(3)

As the ${\Omega }_s$ is constant, $\stackrel{·}{\Phi }={\Omega }_s$ . According to Lagrange’s principle and ignoring the high order terms involving slopes, one can obtain the equations of motion for an asymmetric disk under base angular excitations as follows

$$({\mathbf{m}}_d+{\mathbf{m}}_{\mathrm{dt}}){\stackrel{··}{\mathbf{q}}}_d+({\mathbf{g}}_d+{\mathbf{g}}_{\mathrm{dt}}){\stackrel{·}{\mathbf{q}}}_d+{\mathbf{k}}_{\mathrm{dt}}{\mathbf{q}}_d={\mathbf{f}}_{\mathrm{dg}}+{\mathbf{s}}_{\mathrm{dt}}+{\mathbf{f}}_{\mathrm{du}}$$

(4)

where $m_d,g_d$ is the constant mass and gyroscopic matrices and

$$\begin{array}{l}{m}_d=\mathrm{diag}([m_d\mathrm{,}{m}_d\mathrm{,}{I}_m\mathrm{,}{I}_m]),\\ g_d=\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& I_p{\Omega }_s\\ 0& 0& -I_p{\Omega }_s& 0\end{array}\right]\end{array}$$

(5)

The ${\mathbf{f}}_{\mathrm{dg}}=[\begin{array}{cccc}0\hfill & -m_{d}g\hfill & 0\hfill & 0\hfill \end{array}{]}^T,{\mathbf{f}}_{\mathrm{du}}$ represents the self-gravity force and unbalanced excitations of the disk. In equation (4), the ${\mathbf{m}}_{\mathrm{dt}},{\mathbf{g}}_{\mathrm{dt}},\mathrm{and}{\mathbf{k}}_{\mathrm{dt}}$ are additional mass, gyroscopic, and stiffness matrices induced by the disk asymmetry and time-varying base motions. Moreover, the base motions also induce external excitations to the disk equations, and the expressions are denoted by ${\mathbf{s}}_{\mathrm{dt}}$ . Through simplification and arrangement, one can obtain the expressions for the additional matrices and vectors as follows

$${\mathbf{m}}_{\mathrm{dt}}=\left[\begin{array}{cccc}0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & I_m\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & I_m\hfill \end{array}\right]\left[\begin{array}{cccc}0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & m_{33}\hfill & m_{34}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & m_{43}\hfill & m_{44}\hfill \end{array}\right]$$

(6)

with

$$\begin{array}{c}{m}_{33}=-m_{44}=-{\delta }_d\cos 2{\Omega }_{s}t\hfill \\ m_{34}=m_{43}={\delta }_d\sin 2{\Omega }_{s}t\hfill \end{array}$$

(7)

$${\mathbf{g}}_{\mathrm{dt}}=\left[\begin{array}{cccc}{m}_d\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m_d\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & I_m\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & I_m\hfill \end{array}\right]\left[\begin{array}{cccc}0\hfill & g_{12}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{g_{21}}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & g_{33}\hfill & g_{34}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & g_{43}\hfill & g_{44}\hfill \end{array}\right]$$

(8)

with

$$\begin{array}{c}{g}_{12}=-g_{21}=2{\stackrel{·}{\beta }}_b\hfill \\ g_{33}=-g_{44}=2{\delta }_d{\Omega }_s\sin 2{\Omega }_{s}t\hfill \\ g_{34}=({\chi }_d-2){\stackrel{·}{\beta }}_b+2{\delta }_d{\Omega }_s\cos 2{\Omega }_{s}t\hfill \\ g_{43}=(-{\chi }_d+2){\stackrel{·}{\beta }}_b+2{\delta }_d{\Omega }_s\cos 2{\Omega }_{s}t\hfill \end{array}$$

(9)

$${\mathbf{k}}_{\mathrm{dt}}=\left[\begin{array}{cccc}{m}_d\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m_d\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & I_m\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & I_m\hfill \end{array}\right]\left[\begin{array}{cccc}{k}_{11}\hfill & k_{12}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{k_{21}}& k_{22}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & k_{33}\hfill & k_{34}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & k_{43}\hfill & k_{44}\hfill \end{array}\right]$$

(10)

with

$$\begin{array}{c}k=-{({\stackrel{·}{\beta }}_b^2+{\stackrel{·}{\gamma }}_b^2)}_{11}\hfill \\ k_{12}=({\stackrel{··}{\beta }}_b+{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b)\hfill \\ k_{21}=(-{\stackrel{··}{\beta }}_b+{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b)\hfill \\ k_{22}=-({\stackrel{·}{\alpha }}_b^2+{\stackrel{·}{\beta }}_b^2)\hfill \\ k_{33}={\chi }_d{\Omega }_s{\stackrel{·}{\beta }}_b+2{\delta }_d{\Omega }_s{\stackrel{·}{\beta }}_b\cos 2{\Omega }_{s}t\hfill \\ +{\delta }_d{\stackrel{··}{\beta }}_b\sin 2{\Omega }_{s}t-{\delta }_d{\stackrel{·}{\beta }}_b^2\cos 2{\Omega }_{s}t\hfill \\ +{\stackrel{·}{\beta }}_b^2({\chi }_d-1)+{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b\sin 2{\Omega }_{s}t\hfill \\ +{\delta }_d{\stackrel{·}{\alpha }}_b^2\cos 2{\Omega }_{s}t+(1-{\chi }_d){\stackrel{·}{\alpha }}_b^2\hfill \\ k_{34}=-2{\delta }_d{\Omega }_s{\stackrel{·}{\beta }}_b\sin 2{\Omega }_{s}t+({\chi }_d-1){\stackrel{··}{\beta }}_b\hfill \\ +{\delta }_d{\stackrel{··}{\beta }}_b\cos 2{\Omega }_{s}t+{\delta }_d{\stackrel{·}{\beta }}_b^2\sin 2{\Omega }_{s}t\hfill \\ -{\delta }_d{\stackrel{·}{\alpha }}_b^2\sin 2{\Omega }_{s}t+({\chi }_d-1){\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b\hfill \\ +{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b\cos 2{\Omega }_{s}t\hfill \\ k_{43}=-2{\delta }_d{\Omega }_s{\stackrel{·}{\beta }}_b\sin 2{\Omega }_{s}t+{\delta }_d{\stackrel{··}{\beta }}_b\cos 2{\Omega }_{s}t\hfill \\ +{\stackrel{··}{\beta }}_b+{\delta }_d{\stackrel{·}{\beta }}_b^2\sin 2{\Omega }_{s}t\hfill \\ -{\delta }_d{\stackrel{·}{\alpha }}_b^2\sin 2{\Omega }_{s}t+({\chi }_d-1){\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b+{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b\cos 2{\Omega }_{s}t\hfill \\ k_{44}={\chi }_d{\Omega }_s{\stackrel{·}{\beta }}_b-2{\delta }_d{\Omega }_s{\stackrel{·}{\beta }}_b\cos 2{\Omega }_{s}t\hfill \\ -{\delta }_d{\stackrel{··}{\beta }}_b\sin 2{\Omega }_{s}t+({\chi }_d-1){\stackrel{·}{\beta }}_b^2\hfill \\ +{\delta }_d{\stackrel{·}{\beta }}_b^2\cos 2{\Omega }_{s}t+(1-{\chi }_d){\stackrel{·}{\gamma }}_b^2-{\delta }_d{\stackrel{·}{\gamma }}_b^2\cos 2{\Omega }_{s}t\hfill \\ +{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\gamma }}_b\sin 2{\Omega }_{s}t\hfill \end{array}$$

(11)

$${\mathbf{s}}_{\mathrm{dt}}=\left[\begin{array}{cccc}{m}_d{l}_d\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& m_d{l}_d\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & I_m\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & I_m\hfill \end{array}\right]\left[\begin{array}{c}{s}_1\hfill \\ \multicolumn{1}{c}{s_2}\\ \multicolumn{1}{c}{s_3}\\ \multicolumn{1}{c}{s_4}\end{array}\right]$$

(12)

with

$$\begin{array}{c}{s}_1={\stackrel{··}{\gamma }}_b-{\stackrel{·}{\alpha }}_b{\stackrel{·}{\beta }}_b\hfill \\ s_2=-({\stackrel{··}{\alpha }}_b+{\stackrel{·}{\beta }}_b{\stackrel{·}{\gamma }}_b)\hfill \\ s_3=-{\chi }_d{\Omega }_s{\stackrel{·}{\alpha }}_b-2{\delta }_d{\Omega }_s{\stackrel{·}{\alpha }}_b\cos 2{\Omega }_{s}t\hfill \\ -{\delta }_d{\stackrel{··}{\alpha }}_b\sin 2{\Omega }_{s}t-2{\delta }_d{\Omega }_s{\stackrel{·}{\gamma }}_b\sin 2{\Omega }_{s}t\hfill \\ -{\stackrel{··}{\gamma }}_b+{\delta }_d{\stackrel{··}{\gamma }}_b\cos 2{\Omega }_{s}t-({\chi }_d-1){\stackrel{·}{\alpha }}_b{\stackrel{·}{\beta }}_b\hfill \\ +{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\beta }}_b\cos 2{\Omega }_{s}t\hfill \\ +{\delta }_d{\stackrel{·}{\beta }}_b{\stackrel{·}{\gamma }}_b\sin 2{\Omega }_{s}t\hfill \\ s_4=2{\delta }_d{\Omega }_s{\stackrel{·}{\alpha }}_b\sin 2{\Omega }_{s}t-{\stackrel{··}{\alpha }}_b-{\delta }_d{\stackrel{··}{\alpha }}_b\cos 2{\Omega }_{s}t\hfill \\ +{\chi }_d{\Omega }_s{\stackrel{·}{\gamma }}_b-2{\delta }_d{\Omega }_s{\stackrel{·}{\gamma }}_b\cos 2{\Omega }_{s}t\hfill \\ -{\delta }_d{\stackrel{··}{\gamma }}_b\sin 2{\Omega }_{s}t-{\delta }_d{\stackrel{·}{\alpha }}_b{\stackrel{·}{\beta }}_b\sin 2{\Omega }_{s}t+({\chi }_d-1){\stackrel{·}{\beta }}_b{\stackrel{·}{\gamma }}_b\hfill \\ +{\delta }_d{\stackrel{·}{\beta }}_b{\stackrel{·}{\gamma }}_b\cos 2{\Omega }_{s}t\hfill \end{array}$$

(13)

The unbalanced mass of the disk is denoted by $m_u$ with eccentricity $e$ and phase angle $\phi$ (with respect to the axis $Y$ ). Similar with that of the derivation of disk equations, the kinetic energy of the unbalance mass could be expressed as

$$T_u=\frac{1}{2}{m}_{u}e({\stackrel{·}{\overline{u}}}_u^2+{\stackrel{·}{\overline{v}}}_u^2+{\stackrel{·}{\overline{w}}}_u^2)$$

(14)

where ${\stackrel{·}{\overline{u}}}_u,{\stackrel{·}{\overline{v}}}_u,{\stackrel{·}{\overline{w}}}_u$ are the absolute velocity of the unbalance mass and are given by

$$\begin{array}{c}{\stackrel{·}{\overline{u}}}_u={\stackrel{·}{\beta }}_b(z_d+e\cos ({\Omega }_{s}t+\phi ))-{\stackrel{·}{\gamma }}_b{l}_d+{\stackrel{·}{x}}_d+e{\Omega }_s\cos ({\Omega }_{s}t+\phi )\hfill \\ {\stackrel{·}{\overline{v}}}_u={\stackrel{·}{\gamma }}_b(x_d+e\sin ({\Omega }_{s}t+\phi ))-{\stackrel{·}{\alpha }}_b(z_d+e\cos ({\Omega }_{s}t+\phi ))\hfill \\ {\stackrel{·}{\overline{w}}}_u={\stackrel{·}{\alpha }}_b{l}_d-{\stackrel{·}{\beta }}_b(x_d+e\sin ({\Omega }_{s}t+\phi ))+{\stackrel{·}{z}}_d-e{\Omega }_s\sin ({\Omega }_{s}t+\phi )\hfill \end{array}$$

(15)

According to Lagrange’s principle, the unbalanced mass excitation vector is written as follows

$$f_{du}=m_{u}e\left[\begin{array}{c}\left({({\Omega }_s+{\stackrel{.}{\beta }}_b)}^2+{\stackrel{.}{\gamma }}_b^2\right)\sin ({\Omega }_{s}t+\phi )-({\stackrel{··}{\beta }}_b+{\stackrel{.}{\alpha }}_b{\stackrel{.}{\gamma }}_b)\cos ({\Omega }_{s}t+\phi )\\ \left({({\Omega }_s+{\stackrel{.}{\beta }}_b)}^2+{\stackrel{.}{\alpha }}_b^2\right)\cos ({\Omega }_{s}t+\phi )+({\stackrel{··}{\beta }}_b-{\stackrel{.}{\alpha }}_b{\stackrel{.}{\gamma }}_b)\sin ({\Omega }_{s}t+\phi )\\ 0\\ 0\end{array}\right]$$

(16)

The asymmetric flexible shaft

Two-node Euler beam is used to model the flexible shaft elements. The degrees of freedom include the nodal transverse deflections and slopes in the rotor frame $X-Y-Z$ as ${\mathbf{q}}_s=[\begin{array}{cccccccc}{x}_{s1}\hfill & z_{s1}\hfill & {\theta }_{s1}^z\hfill & {\theta }_{s1}^x\hfill & x_{s2}\hfill & z_{s2}\hfill & {\theta }_{s2}^z\hfill & {\theta }_{s2}^x\hfill \end{array}{]}^T$ . In order to calculate the transverse displacements of a pint inside the shaft element, an element shape function $\mathbf{ѱ}(y)$ should be defined. For the Euler beam element, the expression of $\mathbf{ѱ}(y)$ is given by

$$\mathbf{ѱ}(y)=\left[\begin{array}{cccccccc}{\psi }_1\hfill & 0\hfill & 0\hfill & {\psi }_2\hfill & {\psi }_3\hfill & 0\hfill & 0\hfill & {\psi }_4\hfill \\ \multicolumn{1}{c}{0}& {\psi }_1\hfill & -{\psi }_2\hfill & 0\hfill & 0\hfill & {\psi }_3\hfill & -{\psi }_4\hfill & 0\hfill \end{array}\right]$$

(17)

with

$$\begin{array}{c}{\psi }_1=1-3(y/l{)}^2+2(y/l{)}^3;{\psi }_2=l((y/l)-2(y/l{)}^2+(y/l{)}^3)\hfill \\ {\psi }_3=3(y/l{)}^2-2(y/l{)}^3;{\psi }_4=l(-(y/l{)}^2+(y/l{)}^3)\hfill \end{array}$$

(18)

where $l$ is the element length. The transverse displacements of a pint inside the shaft element are calculated through interpolation as $[x_s(y)z_s(y){]}^T=\mathbf{ѱ}(y){\mathbf{q}}_s$ . Then, one can obtain the total kinetic energy of an asymmetric shaft element as

$$T_s=\underset{0}{\overset{l}{\int }}(\frac{1}{2}\rho A_s({\stackrel{·}{u}}_s^2(y)+{\stackrel{·}{v}}_s^2(y)+{\stackrel{·}{w}}_s^2(y))+\frac{1}{2}\rho (I_{\mathrm{sx}}{\omega }_{\mathrm{sx}}^2+I_{\mathrm{sp}}{\omega }_{\mathrm{sy}}^2+I_{\mathrm{sz}}{\omega }_{\mathrm{sz}}^2))\mathrm{dy}$$

(19)

in which $\rho ,A_s$ are the material density and area of shaft cross section. $I_{\mathrm{sx}}$ and $I_{\mathrm{sz}}$ are the second moments of area about principal axes $x$ and $z$ of the shaft, respectively, and $I_{\mathrm{sp}}$ is the polar moment of area. For flexible shaft with asymmetric cross sections, one has $I_{\mathrm{sx}}\ne I_{\mathrm{sz}}$ . Similar with that of the disk, three parameters are defined: $I_{\mathrm{sm}}=1/2(I_{\mathrm{sx}}+I_{\mathrm{sz}})$ , ${\delta }_s=I_{\mathrm{sx}}-I_{\mathrm{sz}}/2I_{\mathrm{sm}}$ , and ${\chi }_s=I_{\mathrm{sp}}/I_{\mathrm{sm}}$ . The ${\stackrel{·}{u}}_s,{\stackrel{·}{v}}_s,\mathrm{and}{\stackrel{·}{w}}_s$ and ${\omega }_{\mathrm{sx}},{\omega }_{\mathrm{sy}},\mathrm{and}{\omega }_{\mathrm{sz}}$ represent the absolute translational and rotational velocities of the shaft element and could be derived similarly with that of disk (see equations (2) and (3)). By substituting different components of translational and rotational velocities ${\stackrel{·}{u}}_s,{\stackrel{·}{v}}_s,\mathrm{and}{\stackrel{·}{w}}_s$ and ${\omega }_{\mathrm{sx}},{\omega }_{\mathrm{sy}},\mathrm{and}{\omega }_{\mathrm{sz}}$ into equation (19), the detailed expression for the $T_s$ could be obtained. The strain energy due to bending for the asymmetric shaft element is written as

$$\begin{array}{cc}\multicolumn{1}{c}{V_s=}& \frac{1}{2}\underset{0}{\overset{l}{\int }}E{I}_{\mathrm{sm}}[{(x_s^{″})}^2+{(z_s^{″})}^2]\mathrm{dy}\hfill \\ \multicolumn{1}{c}{}& +\frac{1}{2}\underset{0}{\overset{l}{\int }}E{\delta }_s{I}_{\mathrm{sm}}[{(x_s^{″})}^2\cos 2{\Omega }_{s}t-{(z_s^{″})}^2\cos 2{\Omega }_{s}t-2x_s^{″}{z}_s^{″}\sin 2{\Omega }_{s}t]\mathrm{dy}\hfill \end{array}$$

(20)

where $x_s,z_s$ are the displacements of the center of the cross section in fixed coordinates, $E$ is Young’s modulus of the shaft material, and the operator (·)″ represents the second partial derivative of the operand with respect to axis $y$ . According to Lagrange’s principle, the equations of motion for the shaft element considering the expressions of $T_s$ and $V_s$ are obtained as follows

$$\begin{array}{c}({\mathbf{m}}_s+{\mathbf{m}}_{\mathrm{sa}}(t)){\stackrel{··}{\mathbf{q}}}_s+({\mathbf{g}}_s+{\mathbf{g}}_{\mathrm{sa}}(t)+\varsigma {\mathbf{k}}_s){\stackrel{·}{\mathbf{q}}}_s\hfill \\ +({\mathbf{k}}_s+{\mathbf{k}}_{\mathrm{sa}}(t)+{\mathbf{k}}_{\mathrm{sb}}(t)){\mathbf{q}}_s\hfill \\ ={\mathbf{f}}_{\mathrm{sg}}+{\mathbf{s}}_{\mathrm{sb}}(t)\hfill \end{array}$$

(21)

Similar with that of disk model, the constant mass, gyroscopic, and stiffness matrices of the shaft element are denoted by ${\mathbf{m}}_s,{\mathbf{g}}_s,\mathrm{and}{\mathbf{k}}_s$ . The proportional damping coefficient $\varsigma$ is used to reflect the shaft material damping effect. ${\mathbf{f}}_{\mathrm{sg}}=-\rho g{A}_{s}l[\begin{array}{cccccccc}0\hfill & 1/2\hfill & (1/12)l\hfill & 0\hfill & 0\hfill & 1/2\hfill & -(1/12)l\hfill & 0\hfill \end{array}{]}^T$ is the $8\times 1$ gravity force vector. Additional time-dependent mass, gyroscopic, and stiffness matrices of ${\mathbf{m}}_{\mathrm{sa}}(t),{\mathbf{g}}_{\mathrm{sa}}(t),{\mathbf{k}}_{\mathrm{sa}}(t),{\mathbf{k}}_{\mathrm{sb}}(t)$ and external force excitations $s_{\mathrm{sb}}(t)$ are found in the shaft element equations due to the shaft asymmetry and time-variable base motions. Their expressions are given as follows

$$\begin{array}{c}{m}_{\mathrm{sa}}(t)=m_{33}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}-1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}+m_{34}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}\hfill \end{array}$$

(22)

$$\begin{array}{c}{\mathbf{g}}_{\mathrm{sa}}(t)=g_{12}\underset{0}{\overset{l}{\int }}\rho A_s\mathbf{ѱ}{(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{-1}& 0\hfill \end{array}\right]\mathbf{ѱ}(y)\mathrm{dy}+g_{33}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}-1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}\hfill \\ \multicolumn{1}{c}{+g_{34}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}+g_{43}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}}\end{array}$$

(23)

$$\begin{array}{c}{\mathbf{k}}_{\mathrm{sb}}(t)=k_{11}\underset{0}{\overset{l}{\int }}\rho A_s\mathbf{ѱ}{(y)}^T\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\mathbf{ѱ}(y)\mathrm{dy}+k_{12}\underset{0}{\overset{l}{\int }}\rho A_s\mathbf{ѱ}{(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\mathbf{ѱ}(y)\mathrm{dy}\hfill \\ \multicolumn{1}{c}{+k_{21}\underset{0}{\overset{l}{\int }}\rho A_s\mathbf{ѱ}{(y)}^T\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]\mathbf{ѱ}(y)\mathrm{dy}+k_{22}\underset{0}{\overset{l}{\int }}\rho A_s\mathbf{ѱ}{(y)}^T\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\mathbf{ѱ}(y)\mathrm{dy}}\\ \multicolumn{1}{c}{+k_{33}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}+k_{34}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}}\\ \multicolumn{1}{c}{+k_{43}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}+k_{44}\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right]\mathbf{ѱ}\prime (y)\mathrm{dy}}\end{array}$$

(24)

$$\begin{array}{c}{\mathbf{s}}_{\mathrm{sb}}(t)=s_1\underset{0}{\overset{l}{\int }}\rho A_s(y+(i-1)*l)\mathbf{ѱ}{(y)}^T\left[\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]\mathrm{dy}+s_2\underset{0}{\overset{l}{\int }}\rho A_s(y+(i-1)*l)\mathbf{ѱ}{(y)}^T\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\right]\mathrm{dy}\hfill \\ \multicolumn{1}{c}{+s_3\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\right]\mathrm{dy}+s_4\underset{0}{\overset{l}{\int }}\rho I_{\mathrm{sm}}\mathbf{ѱ}\prime {(y)}^T\left[\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]\mathrm{dy}}\end{array}$$

(25)

It is noted that in the above equations, the ${\chi }_d,{\delta }_d$ appeared in the time-variable coefficients $m_{33},m_{34},g_{33},g_{34},g_{43},k_{33},k_{34},k_{43},k_{44},s_3,s_4$ should be substituted by ${\chi }_s,{\delta }_s$ for asymmetric shaft elements. In addition, the time-periodic stiffness matrix ${\mathbf{k}}_{\mathrm{sa}}(t)$ induced by the shaft asymmetry is given by

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathbf{k}}_{\mathrm{sa}}(t)=}& E{\delta }_s{I}_{\mathrm{sm}}\cos 2{\Omega }_{s}t\underset{0}{\overset{l}{\int }}{\mathbf{ѱ}}^{″}{(y)}^T\left[\begin{array}{cc}1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -1\hfill \end{array}\right]{\mathbf{ѱ}}^{″}(y)\mathrm{dy}\hfill \\ \multicolumn{1}{c}{}& +E{\delta }_s{I}_{\mathrm{sm}}\sin 2{\Omega }_{s}t\underset{0}{\overset{l}{\int }}{\mathbf{ѱ}}^{″}{(y)}^T\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{1}& 0\hfill \end{array}\right]{\mathbf{ѱ}}^{″}(y)\mathrm{dy}\hfill \end{array}$$

(26)

Bearing support

Linear flexible bearing model is adopted in this study, and the influence of slopes and bending moment is neglected in the model. The expressions for the stiffness and damping matrices ${\mathbf{k}}_{\mathrm{brg}}$ and ${\mathbf{c}}_{\mathrm{brg}}$ could be found in Nandi and Neogy. ¹⁷

Global FEM equations

Without loss of generality, the harmonic functions are utilized to simulate the time-periodic base motions.^21–23 The angular motions of the base are assumed to vary harmonically with the same frequency ( ${\omega }_b$ )

$$\begin{array}{c}{\alpha }_b={\alpha }_{b0}+{\alpha }_{b1}\cos {\omega }_{b}t\hfill \\ \multicolumn{1}{c}{{\beta }_b={\beta }_{b0}+{\beta }_{b1}\cos {\omega }_{b}t}\\ \multicolumn{1}{c}{{\gamma }_b={\gamma }_{b0}+{\gamma }_{b1}\cos {\omega }_{b}t}\end{array}$$

(27)

where the constant and time-variant amplitudes of the base angular motions are denoted by ${\alpha }_{b0},{\beta }_{b0},{\gamma }_{b0}$ and ${\alpha }_{b1},{\beta }_{b1},{\gamma }_{b1}$ , respectively.

Through assembly of element equations of motion for the asymmetric disk and flexible shaft, the FEM for the whole rotor-bearing system under base angular excitations is obtained as follows

$$(M+M_v(t))\stackrel{··}{\mathbf{q}}+(G_0+C_0+G_v(t))\stackrel{.}{q}+(K_0+K_v(t))q=R(t)+F_u(t)$$

(28)

where $\mathbf{q}=[{\mathbf{q}}_1^T\cdots {\mathbf{q}}_i^T\cdots {\mathbf{q}}_{N_e+1}^T]$ is the $4(N_e+1)\times 1$ ( $N_e$ is the element number) dimension global displacement vector; $\mathbf{M}=\sum ({\mathbf{m}}_d+{\mathbf{m}}_s),{\mathbf{G}}_0=\sum ({\mathbf{g}}_d+{\mathbf{g}}_s),{\mathbf{C}}_0=\sum (\varsigma {\mathbf{k}}_s+{\mathbf{c}}_{\mathrm{brg}}),{\mathbf{K}}_0=\sum ({\mathbf{k}}_s+{\mathbf{k}}_{\mathrm{brg}})$ are the $4(N_e+1)\times 4(N_e+1)$ dimension mass, gyroscopic, damping, and stiffness matrices of the system, and the symbol ∑ represents the assembly of FE matrices after appropriate globalization method is followed. The additional time-varying mass, gyroscopic, and stiffness matrices ${\mathbf{M}}_v(t),{\mathbf{G}}_v(t),{\mathbf{K}}_v(t)$ and external force vectors $R(t),F_u(t)$ are expressed as

$$\begin{array}{c}{\mathbf{M}}_v(t)=\sum ({\mathbf{m}}_{\mathrm{dt}}+{\mathbf{m}}_{\mathrm{sa}}(t))\hfill \\ {\mathbf{G}}_v(t)=\sum ({\mathbf{g}}_{\mathrm{dt}}+{\mathbf{g}}_{\mathrm{sa}}(t))\hfill \\ {\mathbf{K}}_v(t)=\sum ({\mathbf{k}}_{\mathrm{dt}}+{\mathbf{k}}_{\mathrm{sa}}(t)+{\mathbf{k}}_{\mathrm{sb}}(t))\hfill \\ \mathbf{R}(t)=\sum ({\mathbf{s}}_{\mathrm{dt}}+{\mathbf{s}}_{\mathrm{sb}}(t))\hfill \\ {\mathbf{F}}_u(t)=\sum ({\mathbf{f}}_{\mathrm{du}}+{\mathbf{f}}_{\mathrm{dg}}+{\mathbf{f}}_{\mathrm{sg}})\hfill \end{array}$$

(29)

Without base angular motions, the ${\mathbf{M}}_v(t),{\mathbf{G}}_v(t),{\mathbf{K}}_v(t)$ of asymmetric rotor system are all time-periodic with the frequency of $2{\Omega }_s$ . After considering the harmonic base motions, besides the $2{\Omega }_s$ , the time-varying elements of the above matrices contain multiple frequencies, including $2{\Omega }_s\pm {\omega }_b,2{\Omega }_s\pm 2{\omega }_b,{\omega }_b,2{\omega }_b$ . This can be proved from equations (9) and (11) through simple trigonometric calculations. The external force vectors $\mathbf{R}(t),{\mathbf{F}}_u(t)$ also contain multiple frequencies as long as the base angular motions are taken into account.

Computations and discussions

Model validation

The rotor model of Oncescu et al. ¹⁶ is adopted here. The disk is symmetric ( ${\delta }_d=0$ ) and the shaft is asymmetric ( ${\delta }_s\ne 0$ ). The parameters of the rotor system are $m_d=2\mathrm{kg}$ , $I_m=0.005\mathrm{kg}{\mathrm{m}}^2$ , ${\chi }_d=2$ , rotor length $L=1\mathrm{m}$ , $\rho =7750\mathrm{kg}{\mathrm{m}}^3$ , $E=2\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$ , $I_{\mathrm{sm}}=4\times {10}^{-8}{\mathrm{m}}^4$ , ${\chi }_s=2$ , and $A_s=0.693\times {10}^{-3}{\mathrm{m}}^2$ . The shaft is divided into five elements and the disk locates at node 3. The proportional damping of the shaft is neglected, that is, $\varsigma =0$ . The stiffness and damping coefficients of left bearing are $k_{\mathrm{xx}}=3.5\times {10}^5\mathrm{N}/\mathrm{m}$ , $k_{\mathrm{zz}}=5.5\times {10}^5\mathrm{N}/\mathrm{m}$ , $k_{\mathrm{xz}}=-1\times {10}^5\mathrm{N}/\mathrm{m}$ , $c_{\mathrm{xx}}=26\mathrm{N}\mathrm{s}/\mathrm{m}$ , $c_{\mathrm{zz}}=20\mathrm{N}\mathrm{s}/\mathrm{m}$ , and $c_{\mathrm{xz}}=-8\mathrm{N}\mathrm{s}/\mathrm{m}$ . Similarly, right bearing’s parameters are $k_{\mathrm{xx}}=4.5\times {10}^5\mathrm{N}/\mathrm{m}$ , $k_{\mathrm{zz}}=3\times {10}^5\mathrm{N}/\mathrm{m}$ , $k_{\mathrm{xz}}=-1\times {10}^5\mathrm{N}/\mathrm{m}$ , $c_{\mathrm{xx}}=24\mathrm{N}\mathrm{s}/\mathrm{m}$ , $c_{\mathrm{zz}}=30\mathrm{N}\mathrm{s}/\mathrm{m}$ , and $c_{\mathrm{xz}}=-4\mathrm{N}\mathrm{s}/\mathrm{m}$ .

Without considering the base angular motions, the parametric instability induced by the asymmetric rotor could be computed numerically utilizing the discrete state transition matrix (DSTM) method. ³⁴ Both the obtained instability regions and the results in Oncescu et al. ¹⁶ are shown in Figure 2. It is shown that the obtained instability regions of this study are in good agreement with those of Oncescu et al., ¹⁶ indicating that the asymmetric rotor model derived in this study is believable.

Figure 2.

Instability comparisons for the asymmetric rotor-bearing system without base angular motions: (a) results of Oncescu et al. ¹⁶ and (b) results of this study.

The dynamic response of the symmetric rotor-bearing system ( ${\delta }_d={\delta }_s=0$ ) obtained by numerical integration method (Runge–Kutta method) is verified with the results of Driot et al. ²¹ The shaft is divided into four elements, and the disk locates at node 4. The bearing stiffness coefficients are $k_{\mathrm{xx}}=k_{\mathrm{zz}}=1\times {10}^9\mathrm{N}/\mathrm{m}$ to simulate the simply supported boundary condition. The rotor parameters are set to be the same with Driot et al. ²¹ The pitching base motion with ${\alpha }_{b1}=2\times {10}^{-5}\mathrm{rad}$ and three values of frequency ${\omega }_b=251.2,376.8,\mathrm{and}502.4\mathrm{rad}/\mathrm{s}(40,60,\mathrm{and}80\mathrm{Hz})$ are, respectively, considered in the analysis. The rotating speed is ${\Omega }_s=125.6\mathrm{rad}/\mathrm{s}(20\mathrm{Hz})$ and the shaft center orbits of the disk nodes are given in Figure 3. The corresponding results of Driot et al. ²¹ are also presented in Figure 3. It is shown that the ranges of the trajectory obtained by this study are consistent with the reference ones, while slight differences could also be found in the shapes of trajectory. In Driot et al., ²¹ the sinusoidal shape function is used to obtain the 2-degree-of-freedom rotor model. As the disk does not locate at the mid-span, such treatment might cause the differences.

Figure 3.

Shaft center orbit comparisons for the rotor-bearing system under pitching base motion $({\alpha }_{b1}=2\times {10}^{-5})$ and ${\Omega }_s=125.6\mathrm{rad}/\mathrm{s}$ : (a, c, and e) results of Driot et al. ²¹ and (b, d, and f) results of this study: (a) ω_b = 251.2 rad/s, (b) ω_b = 251.2 rad/s, (c) ω_b = 376.8 rad/s, (d) ω_b = 376.8 rad/s, (e) ω_b = 502.4 rad/s, and (f) ω_b = 502.4 rad/s.

Parametric instability

In the following analysis, the rotor model of Oncescu et al. ¹⁶ is adopted, and the only difference is that the disk location is moved to node 2. Before instability analysis, the first order of whirling frequency varying with rotating speed is shown in Figure 4. Both the asymmetry and base angular motions are neglected. The first order of backward and forward critical frequencies are ${\omega }_{c1}=220.4\mathrm{rad}/\mathrm{s}$ and ${\omega }_{c2}=251.3\mathrm{rad}/\mathrm{s}$ , which are also marked in Figure 4.

Figure 4.

Whirling analysis of the symmetric rotor-bearing system without base angular motions.

Without base movements, instability regions of the asymmetric rotor system are obtained numerically by the DSTM method and presented in Figure 5. Three values of proportional damping of the shaft, that is, $\varsigma =0,1\times {10}^{-5},5\times {10}^{-5}\mathrm{s}$ , are considered in the analysis. It is shown that three instability regions, of which the starting points are about ${\omega }_{c1},{\omega }_{c2},({\omega }_{c1}+{\omega }_{c2})/2$ , appear in Figure 5. They are identified as the two primary instability regions and one combination instability region according to the parametric vibration theory. ³⁵ As long as the damping is considered, all the instability regions are reduced gradually. Especially for the lower asymmetries $({\delta }_s<0.05)$ , the instability regions have been vanished.

Figure 5.

Instability regions for the asymmetric rotor system without base movements: (a) $\varsigma =0$ , (b) $\varsigma =1\times {10}^{-5}$ , and (c) $\varsigma =5\times {10}^{-5}$ .

Dynamic response

The forced steady-state response of the system under unbalanced excitations, shaft asymmetry, and base angular motions will be analyzed in this section. If the system is in unstable region, the response amplitudes will be enlarged exponentially with time under initial conditions. Obviously, in this case, study on the dynamic response of the system does not make any sense. Thus, the system should be stable in dynamic response analysis. According to the instability regions in Figure 5, the shaft asymmetry and damping are chosen as ${\delta }_s=0.05$ and $\varsigma =5\times {10}^{-5}\mathrm{s}$ . The unbalanced mass excitations are $m_{u}e=4\times {10}^{-3}\mathrm{kg}\mathrm{m}$ and $\phi =\pi /4$ . In the following, the dynamic response of the system under pitching, rolling, and yawing motions are analyzed.

Pitching base motion is considered

Without base movements, the contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system is presented in Figure 6. One can see that the response spectra mainly contain the rotating speed and its integral multiples $({\Omega }_s,2{\Omega }_s)$ and natural frequencies $({\omega }_{c1},{\omega }_{c2})$ . It is known that the spectral amplitudes of ${\omega }_{c1},{\omega }_{c2}$ would be attenuated gradually with time in damped systems except for the resonant system. Besides the backward and forward critical resonances ${\omega }_{c1},{\omega }_{c2}$ , the rotating speed equal to half of the critical speeds ${\omega }_{c1}/2,{\omega }_{c2}/2$ would also lead to the external resonances, as shown in Figure 6. This is due to the parametric stiffness excitation from shaft asymmetry. ¹⁶

Figure 6.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system without base movements.

When the periodic pitching motion is considered ( ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ and ${\alpha }_{b1}=0.05\mathrm{and}0.1\mathrm{rad}$ ), the contour plots of the steady-state response are given in Figures 7 and 8. Both the response spectra and external resonances are greatly changed after the pitching motion is considered. For the response spectra, besides ${\Omega }_s,2{\Omega }_s$ and ${\omega }_{c1},{\omega }_{c2}$ , there are also base frequency ${\omega }_b$ and combined frequencies between twice of rotating speed and base frequency, that is, $2{\Omega }_s\pm {\omega }_b$ . Increasing the pitching amplitudes ${\alpha }_{b1}(0.05\to 0.1\mathrm{rad})$ makes these spectra more obvious, as shown in Figure 8. As long as the pitching motion is in operation, the spectra of natural frequencies ${\omega }_{c1},{\omega }_{c2}$ become significant for the entire rotation speed range (see Figures 7 and 8). This indicates that the operation of pitching motion causes the transient response more severe and less likely to be weakened by damping.

Figure 7.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic pitching motion ( ${\alpha }_{b1}=0.05\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

Figure 8.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic pitching motion ( ${\alpha }_{b1}=0.1\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

Besides the critical and half-critical resonances, there are also multiple resonances around the half-critical resonances with the interval of half of base frequency ${\omega }_b/2$ . They could be expressed by $({\omega }_{c1}\pm {\omega }_b)/2$ and $({\omega }_{c2}\pm {\omega }_b)/2$ , as shown in Figures 7 and 8. This is the result of combined effects of the periodic pitching motion and shaft asymmetry. Also, increasing the value of ${\alpha }_{b1}$ makes the external resonances more distinct.

Rolling base motion is considered

For the operation of the periodic rolling motion ( ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ and ${\beta }_{b1}=0.1\mathrm{and}0.2\mathrm{rad}$ ), Figures 9 and 10 show the contour plots of the steady-state response. One can also find multiple response spectra in the dynamic response after the rolling motion is considered. Different with that of the pitching motion, the additional response spectra are around the rotation speed, that is, ${\Omega }_s\pm {\omega }_b$ . Moreover, the spectrum of base frequency ${\omega }_b$ is not found and the amplitudes of natural frequencies ${\omega }_{c1},{\omega }_{c2}$ are not as obvious as that of the pitching base–excited system. This indicates that the transient response would not be significantly excited by the periodic rolling base motion. Increasing the base amplitude ${\beta }_{b1}$ from 0.1 rad to 0.2 rad not only causes the current spectra more distinct, but also induces more spectra around the rotation speed and twice of rotation speed. They are identified by ${\Omega }_s\pm 2{\omega }_b$ and $2{\Omega }_s\pm {\omega }_b$ , as shown in Figure 10.

Figure 9.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic rolling motion ( ${\beta }_{b1}=0.1\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

Figure 10.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic rolling motion ( ${\beta }_{b1}=0.2\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

The periodic rolling motion also induces multiple external resonances besides the critical and half-critical resonances. However, different from that of the pitching base–excited system, the additional resonances of rolling base–excited system are around the critical resonances with the interval of base frequency ${\omega }_b$ . They are identified by ${\omega }_{c1}\pm {\omega }_b$ and ${\omega }_{c2}\pm {\omega }_b$ . It is noted that the resonance of ${\omega }_{c2}+{\omega }_b$ is greater than 300 rad/s, and thus, it is not marked in the Figures 9 and 10. Moreover, the increase in base amplitudes ${\beta }_{b1}$ more obviously leads to the external resonances.

Yawing base motion is considered

When the periodic yawing motion is considered ( ${\gamma }_{b1}=0.05,0.1\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ), the results are shown in Figures 11 and 12. One can find that both the response spectra and external resonant characteristics are similar with that of pitching base–excited system. The response spectra are summarized as ${\Omega }_s,2{\Omega }_s$ , ${\omega }_{c1},{\omega }_{c2}$ (which are also found in asymmetric system), and ${\omega }_b,2{\Omega }_s\pm {\omega }_b$ . The multiple external resonances are expressed as ${\omega }_{c1},{\omega }_{c2}$ , ${\omega }_{c1}/2,{\omega }_{c2}/2$ (which are also found in asymmetric system), and $({\omega }_{c1}\pm {\omega }_b)/2,({\omega }_{c2}\pm {\omega }_b)/2$ . Similarly, increasing the base amplitude ${\gamma }_{b1}$ causes both the response spectra and additional resonances more distinct.

Figure 11.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic yawing motion ( ${\gamma }_{b1}=0.05\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

Figure 12.

Contour plot of the steady-state response (x-direction of node 3) for the asymmetric rotor system under periodic yawing motion ( ${\gamma }_{b1}=0.1\mathrm{rad}$ and ${\omega }_b=50\mathrm{rad}/\mathrm{s}$ ).

Conclusion

Dynamic response characteristics of asymmetric rotor systems under time-periodic base angular motions are studied numerically in this article. As long as the periodic base angular motion is considered, both the response spectra and external resonances change significantly.

For the pitching or yawing motion in operation, additional frequencies, including base frequency and combined frequencies between twice of rotating speed and base frequency, are found in the dynamic response. The spectra of natural frequencies become significant for the entire rotation speed range, indicating that the operation of pitching or yawing motion causes the transient response more severe and less likely to be weakened by damping. When the periodic rolling motion is considered, the additional spectra observed in the dynamic response are the combined frequencies between the rotation speed and base frequency, which are different from that of the pitching or yawing base-excited system. In addition, the base frequency does not exist in the response, and natural frequencies are not obvious.

Besides the critical and half-critical resonances (which are found in the asymmetric rotor system), there are multiple resonances after considering the periodic base motions. For the pitching or yawing base–excited system, the multiple resonances are around the half-critical resonances with the interval of half of base frequency. However, when the rolling motion is in operation, the multiple resonances are found in the vicinity of critical resonances, and the interval is base frequency.