JZUS - Journal of Zhejiang University SCIENCE

Journal of Zhejiang University SCIENCE A 2014 Vol.15 No.1 P.39-52

Numerical analysis of a nonlinear double disc rotor-seal system*

Author(s): Wen-jie Zhou¹, Xue-song Wei¹, Xian-zhu Wei², Le-qin Wang¹
Affiliation(s): 1. Institute of Process Equipment, Zhejiang University, Hangzhou 310027, China; more
Corresponding email(s): zhouwenjiezwj@zju.edu.cn
Key Words: Nonlinear, Rotor-seal system, Finite element method (FEM), Fluid excitation

Share this article to： More <<< Previous Article \|Next Article >>>

Wen-jie Zhou, Xue-song Wei, Xian-zhu Wei, Le-qin Wang. Numerical analysis of a nonlinear double disc rotor-seal system[J]. Journal of Zhejiang University Science A, 2014, 15(1): 39-52.

@article{title="Numerical analysis of a nonlinear double disc rotor-seal system",
author="Wen-jie Zhou, Xue-song Wei, Xian-zhu Wei, Le-qin Wang",
journal="Journal of Zhejiang University Science A",
volume="15",
number="1",
pages="39-52",
year="2014",
publisher="Zhejiang University Press & Springer",
doi="10.1631/jzus.A1300230"
}

%0 Journal Article
%T Numerical analysis of a nonlinear double disc rotor-seal system
%A Wen-jie Zhou
%A Xue-song Wei
%A Xian-zhu Wei
%A Le-qin Wang
%J Journal of Zhejiang University SCIENCE A
%V 15
%N 1
%P 39-52
%@ 1673-565X
%D 2014
%I Zhejiang University Press & Springer
%DOI 10.1631/jzus.A1300230

TY - JOUR
T1 - Numerical analysis of a nonlinear double disc rotor-seal system
A1 - Wen-jie Zhou
A1 - Xue-song Wei
A1 - Xian-zhu Wei
A1 - Le-qin Wang
J0 - Journal of Zhejiang University Science A
VL - 15
IS - 1
SP - 39
EP - 52
%@ 1673-565X
Y1 - 2014
PB - Zhejiang University Press & Springer
ER -
DOI - 10.1631/jzus.A1300230

Abstract
Chinese Summary
Academic Network
Reviewer Comment

Abstract: Based on the finite element method (FEM) and the Lagrange equation, a novel nonlinear model of a double disc rotor-seal system, including the coupled effects of the gravity force of the discs, Muszynska’s nonlinear seal fluid dynamic force, and the mass eccentricity of the discs, is proposed. The fourth order Runge-Kutta method is applied to solve the motion equations of the system and numerically determine the vibration response of the center of the discs. The dynamic behavior of the system is analyzed using bifurcation diagrams, time-history diagrams, axis orbit diagrams, Poincaré maps, and amplitude spectrums. With the rotor speed increasing, the system presents rich forms including periodic, multi-periodic, quasi-periodic, and chaotic motion. We also discuss the effects of the distance between the two discs, the mass of the discs, seal clearance, seal length, and seal drop pressure on the dynamic behavior of the system. The numerical results demonstrate that a symmetrical disc structure, small disc mass, proper seal clearance, long seal length and high seal drop pressure can enhance the stability of a double disc rotor-seal system. The results provide a theoretical foundation for the design of multi-stage sealing systems.

非线性双圆盘转子-密封系统的数值分析

研究目的：求解双圆盘转子-密封系统的非线性振动特性和运动响应
创新方法：采用有限元法（FEM）和拉格朗日方程求解双圆盘转子-密封系统，进而为研究多级转子系统的非线性振动问题提供有效方法。
研究手段：基于有限元法（FEM）和拉格朗日方程，得到包含Muszynska非线性密封流体力和圆盘重力作用下的双圆盘转子-密封系统的运动方程。同时利用四阶龙格-库塔法求解系统动特性运动响应情况，利用分岔图、时间历程图、轴心轨迹图、庞加莱映射和幅值谱等分析图研究双圆盘转子-密封系统的非线性振动特性。
重要结论：随着转速的增大，双圆盘转子-密封系统呈现丰富的非线性运动形式，包括周期性运动、多周期运动、准周期运动以及混沌运动。在右端圆盘不平衡质量小于34 kg、密封间隙范围为0.376 mm–0.54 mm、密封长度大于0.13 m或者密封压差高于0.104 MPa的情况下均有利于提高双圆盘转子-密封系统的稳定性。

关键词：非线性；转子-密封系统；有限元法；流体激励

Darkslateblue:Affiliate; Royal Blue:Author; Turquoise:Article

Article Content

1. Introduction

External excitations, such as the oil-film force of journal bearings or the fluid exciting force inside the seal structure between a stator and rotor, are important factors that may induce instability in modern rotating machinery. The effect of the fluid exciting force in a multi-stage centrifugal pump, shrouded turbine, steam turbine, and other rotating machinery that includes a multi-level sealing structure can be very significant (Yuan et al., ; Megerle et al., ). Research on the mechanisms of fluid-solid interaction and the influence of nonlinear seal force on rotor-seal systems has become a priority for the safe operation of rotating machinery.

Many experimental studies have been carried out on the dynamic characteristics of various seals (Kaneko et al., ; Smalley et al., ; Childs et al., ). Muszynska () and Muszynska and Bently () proposed a nonlinear seal fluid dynamic force model that considered the circulating velocity as the key factor affecting the stability of the rotor system. Fei et al. () developed a procedure, based on the finite element method (FEM), which can calculate the dynamics of dual rotor systems and obtain the coupling motion equations of the subsystems. Li et al. () applied the Hamilton principle and the FEM to establish a new dynamic model of a rotor system. The model was used to analyze the dynamic behavior of a rotor system based on the Timoshenko model, with the coupled effects of the nonlinear oil film force, the nonlinear seal force, and the mass eccentricity of the disc. Ding et al. () investigated the Hopf bifurcation behavior of a symmetric rotor-seal system using Muszynska’s nonlinear seal fluid dynamic force model. They discovered that for a perfectly balanced system, the instability from certain critical equilibrium positions proved to be the result of Hopf bifurcation and only the supercritical type was found for a specific rotor system using Poore’s algebraic criteria. Li et al. () analyzed the influence of a labyrinth seal on the nonlinear characteristics of a rotor-labyrinth seal system with a sliding bearing, and judged the system’s stability using the Floquet theory and the bifurcation theorem. Wang and Wang () researched the nonlinear coupling vibration and bifurcation of a high-speed centrifugal compressor with a labyrinth seal and two air-film journal bearings under different conditions of pressure drop and seal length. The Muszynska nonlinear seal model was also used successfully by Hua et al. () to investigate the Hopf bifurcation and stability of a rotor-seal system using a high-precision direct integration method. Wang et al. (; ) established a nonlinear mathematical model for orbital motion of a rotor under the influence of leakage flow through a labyrinth seal. They found that the destabilization speed of the rotor was reduced due to the aerodynamic force induced by the leakage flow through the interlocking seal.

Although previous studies paid much attention to the dynamic behavior of rotor-seal or rotor-bearing-seal systems, only the Jeffcott rotor system was considered. No model was developed for a nonlinear double disc rotor-seal system with a coupled effect between the two discs based on the Lagrange equation and FEM. There is an urgent need to research the mechanisms of interaction of different levels of the impeller seal force in multi-stage sealing systems. In this paper, a nonlinear double disc rotor-seal system model is established, based on the Lagrange equation and FEM. Equations of motion, considering the Muszynska nonlinear seal force, are solved by the fourth-order Runge-Kutta method. The nonlinear characteristics of the double disc rotor-seal system are presented in the form of bifurcation diagrams, time-history diagrams, axis orbit diagrams, Poincaré maps, and amplitude spectrums. The effects of the distance between the two discs, the mass of the discs, seal clearance, seal length and seal drop pressure on the dynamic behavior of the system are also studied, and linear and nonlinear seal models are compared.

2. Model of a double disc rotor-seal system

A typical double disc rotor-seal system, consisting of a rotating shaft, wheel, nonlinear seal force and clamped support at each end of the shaft, was modeled (Fig. 1), where L ₁ is the length between the left edge of the shaft and the left disc, L ₂ is the length between the left and right discs, L ₃ is the length between the right disc and the right edge of the shaft, d _d is the diameter of the disc, m _d1 is the mass of the left disc, and m _d2 the mass of the right disc.

Fig.1
Double disc rotor-seal system

A finite element model of the system was established (Fig. 2). The system is divided into a total of four nodes and three shaft elements. Each node has four degrees of freedom, including two translational and two rotational motions. Nodes 1 and 4 are the supporting points, which are allowed only rotational motions, and nodes 2 and 3 bear both the gravity force of the disc and the nonlinear seal force, where F _sx is the seal force in the x-direction, and F _sy the seal force in the y-direction. The subscripts 1 and 2 denote the left and right discs, respectively, where the subscript ‘s’ denotes the seal.

Fig.2
Finite element model of a double disc rotor-seal system

3. Component equations and system motion equations

Considering the double disc rotor-seal system is a multi-degree of freedom and non-conservative system, the Lagrange equation was applied to establish the motion equation of the double disc rotor-seal system. The non-conservative system Lagrange equation usually can be described in the following form: (1) $\frac{{\text{d}}}{{{\text{d}}t}}\left( {\frac{{\partial T}}{{\partial {{\dot u}_j}}}} \right) - \frac{{\partial (T - V)}}{{\partial {u_j}}} = {Q_j}$ , where T is the kinetic, V is the strain energy, u_j are independent generalized coordinates, and Q_j are generalized forces.

3.1. Rotating coordinate

The kinetic energy and strain energy of the components, shaft units and disc, can be derived using the projected angle method (Li et al., ). Fig. 3 shows a schematic map of a rotating coordinate and A-xyz is the fixed coordinate. Because of the unbalanced force of the disc, at a certain rotation speed Ω, deformation of the shaft would occur, causing the center line of the disc to be no longer collinear with the undeformed rotor center line. A rotating coordinate system o′-ξηζ, is attached to the disc and rotates synchronously with the rotor. θ is the angle between the ζ and z axes, and θ_x , θ_y are the angular displacements, defined as the angles between the z axis and the projection of the ζ axis onto the x-z and y-z planes, respectively. L is the length of the shaft, L _d is the place of disc center in the z direction, and r is the distance of disc center in x-y plane.

Fig.3
Schematic map of rotating coordinates

The absolute angular speed ω in the o′-ξηζ coordinate can be transformed to be rotation speed, expressed by the A-xyz coordinate (Zhong, ). The relationship between the fixed A-xyz coordinate and the o′-ξηζ coordinate attached to the disc can be described as follows: (2) ${\kern 0pt} \left\{ \begin{matrix} {\omega _\xi } \hfill \\ {\omega _\eta } \hfill \\ {\omega _\zeta } \hfill \\ \end{matrix} \right\} = \left\{ \begin{matrix} {{\dot \theta }_\xi }\cos \phi + {{\dot \theta }_y}\cos {\theta _\xi }\sin \phi \hfill \\ - {{\dot \theta }_\xi }\sin \phi + {{\dot \theta }_y}\cos {\theta _\xi }\cos \phi \hfill \\ {\kern 29pt} \dot \phi - {{\dot \theta }_y}\sin {\theta _\xi } \hfill \\ \end{matrix} \right\}$ , where θ_ξ , ϕ are relative rotation angles between the o′-ξηζ coordinate and the A-xyz coordinate when the ξ-axis and ζ-axis are respectively, fixed in turn. ω_ξ , ω_η , and ω_ζ are the components of ω in o′-ξηζ coordinate.

Considering θ_ξ ≈θ_x and $\dot \phi \approx \Omega$ , Eq. (2) can be written as (3) ${\kern 0pt} \left\{ \begin{matrix} {\omega _\xi } \hfill \\ {\omega _\eta } \hfill \\ {\omega _\zeta } \hfill \\ \end{matrix} \right\} = \left\{ \begin{matrix} {\kern 4pt} \cos \phi {\kern 9pt} \cos {\theta _x}\sin \phi {\kern 9pt} 0 \hfill \\ - \sin \phi {\kern 7pt} \cos {\theta _x}\cos \phi {\kern 9pt} 0 \hfill \\ {\kern 12pt} 0{\kern 23pt} - \sin {\theta _x}{\kern 12pt} \;{\kern 1pt} 1 \hfill \\ \end{matrix} \right\}\left\{ \begin{matrix} {{\dot \theta }_x} \hfill \\ {{\dot \theta }_y} \hfill \\ \Omega \hfill \\ \end{matrix} \right\}$ .

3.2. Rigid disc

For a rigid disc, the strain energy would be ignored and its degrees of freedom can be expressed by two translational displacements x, y and two rotational displacements θ_x , θ_y . Let ${u_{dx}} = {[x,{\theta _y}]^{\text{T}}},\;\;{u_{dy}} = {[y, - {\theta _x}]^{\text{T}}}$ , then the kinetic energy of the disc for lateral motion is given by (Zhong, ) (4) ${T_d} = \frac{1}{2}\dot u_{dx}^T{M_d}{\dot u_{dx}} + \frac{1}{2}\dot u_{dy}^{\text{T}}{M_d}{\dot u_{dy}}{\kern 1pt} + \Omega \dot u_{dx}^TJ{\dot u_{dy}} + \frac{1}{2}{J_p}{\Omega ^2}$ , where ${\kern 0pt} {M_d} = \left[ \begin{matrix} {m_d}{\kern 9pt} 0 \hfill \\ {\kern 2pt} 0{\kern 10pt} {J_d} \hfill \\ \end{matrix} \right]$ , ${\kern 0pt} J = \left[ \begin{matrix} 0{\kern 12pt} 0 \hfill \\ 0{\kern 8pt} {J_p} \hfill \\ \end{matrix} \right]$ , where J _d is the diametric moment, and J _p is the polar moment of inertia of the disc.

Substituting Eq. (4) into Eq. (1), the motion equations of a rigid disc can be obtained as (5)

${\kern 0pt} \left\{ \begin{matrix} {M_d}{{\dot u}_{dx}} + \Omega J{{\dot u}_{dy}} = {q_{dx}}, \hfill \\ {M_d}{{\dot u}_{dy}} - \Omega J{{\dot u}_{dx}} = {q_{dy}}, \hfill \\ \end{matrix} \right.$

where the subscript ‘d’ denotes the disc and q _dx, q _dy are generalized forces corresponding with u _dx, u _dy, respectively.

3.3. Elastic shaft element

Fig. 4 shows an elastic shaft element, including two nodes A and B, i.e., with eight degrees of freedom. Define the generalized coordinates ${u_{ex}} = {[{x_A},{\theta _{yA}},{x_B},{\theta _{yB}}]^{\text{T}}}$ , and ${u_{ey}} = {[{y_A}, - {\theta _{xA}},{y_B}, - {\theta _{xB}}]^{\text{T}}}$ .

Fig.4
Elastic shaft element

The generalized coordinates of a typical point internal to the element can be described by the endpoint of the shaft (Nelson and McVaugh, ; Zhang, ): (6) $x = N{u_{ex}},\;\;y = N{u_{ey}}$ , (7) ${\theta _y} = N'{u_{ex}},\;\; - {\theta _x} = N'{u_{ey}}$ , where $N = [{N_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {N_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {N_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {N_4}]$ , $N' = \left[ {{{N'}_1}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{N'}_2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{N'}_3}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{N'}_4}} \right]$ , N ₁=1−3(s/l)²+2(s/l)³, N ₂=l[s/l−2(s/l)²+(s/l)³], N ₃=3(s/l)²−2(s/l)³], and N ₄=l[−(s/l)²+(s/l)³].

Similarly, the kinetic energy and strain energy also can be described as a function of the displacements and velocity of the elastic shaft element nodes and obtained by computing an integral over the length of the element (Zhong, ). The kinetic energy and strain energy can be calculated as follows: (8) ${T_e} = \frac{1}{2}\dot u_{ex}^{\text{T}}({M_{eT}} + {M_{eR}}){\dot u_{ex}} + \frac{1}{2}\dot u_{ey}^{\text{T}}({M_{eT}} + {M_{eR}}){\dot u_{ey}}{\kern 1pt} + \Omega \dot u_{ex}^{\text{T}}{J_e}{\dot u_{ey}} + \frac{1}{2}{j_{pe}} \cdot l \cdot \Omega$ , (9) ${V_e} = \frac{1}{2}u_{ex}^{\text{T}}{K_e}{u_{ex}} + \frac{1}{2}u_{ey}^{\text{T}}{K_e}{u_{ey}}$ , where ${M_{eT}} = \int_0^l {\mu {N^{\text{T}}}N} {\text{d}}z,\;\;\;{J_e} = \int_0^l {{j_p}{{N'}^{\text{T}}}N'} {\text{d}}z,$ ${M_{eR}} = \int_0^l {{j_d}{{N'}^{\text{T}}}N'} {\text{d}}z,\;\;\;{K_e} = \int_0^l {EI{{N''}^{\text{T}}}N''} {\text{d}}z,$ where μ is the mass of shaft element unit length.

Substituting Eqs. (8) and (9) into Eq. (1), the motion equations of the elastic shaft element can be obtained as (10)

${\kern 0pt} \left\{ \begin{matrix} {M_e}{{\ddot u}_{ex}} + \Omega {J_e}{{\dot u}_{ey}} + {K_e}{u_{ex}} = {q_{ex}}, \hfill \\ {M_e}{{\ddot u}_{ey}} - \Omega {J_e}{{\dot u}_{ex}} + {K_e}{u_{ey}} = {q_{ey}}, \hfill \\ \end{matrix} \right.$

where subscript ‘e’ denotes the elastic shaft element and q _ex, q _ey are generalized forces corresponding with u _ex, u _ey, respectively, which are always ignored.

3.4. Nonlinear seal force model

Muszynska’s nonlinear seal fluid dynamic force model was adopted to describe the nonlinear seal force acting on the discs. Muszynska’s model has been widely used to describe the nonlinear characteristic of seal force and can be expressed as follows (Muszynska, ; Hua et al., ): (11) ${\kern 0pt} \left\{ \begin{matrix} {F_{sx}} \hfill \\ {F_{sy}} \hfill \\ \end{matrix} \right\} = - \left[ \begin{matrix} {K_s} - {m_s}{\gamma ^2}{\Omega ^2}\;\;\;\;\;\;\;\gamma \Omega {C_s} \hfill \\ \;\;\;\; - \gamma \Omega {C_s}\;\;\;\;\;{K_s} - {m_s}{\gamma ^2}{\Omega ^2}\; \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} x \hfill \\ y \hfill \\ \end{matrix} \right\} - \left[ \begin{matrix} \;\;\;\;\;\;{C_s}\;\;\;\;\;\;\;2{m_s}\gamma \Omega \hfill \\ - 2{m_s}\gamma \Omega \;\;\;\;\;\;\;\;{C_s}\; \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\dot x} \hfill \\ {\dot y} \hfill \\ \end{matrix} \right\} - \left[ \begin{matrix} {m_s}\;\;\;\;0 \hfill \\ \;0\;\;\;\;{m_s} \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\ddot x} \hfill \\ {\ddot y} \hfill \\ \end{matrix} \right\}$ , where ${K_s} = {K_0}{(1 - {\varepsilon ^2})^{ - {n_1}}}$ , ${C_s} = {C_0}{(1 - {\varepsilon ^2})^{ - {n_1}}}$ , $\gamma = {\gamma _0}{(1 - \varepsilon )^{{n_2}}}$ , n ₁=0.5–3, n ₂=0–1, γ ₀<0.5, ε=(x ²+y ²)^1/2/c, c is the radial clearance of the seal, and K ₀, C ₀, and m _s are expressed through the short bearing model of Childs ().

3.5. System motion equations

For a double disc rotor-seal system, the motion equations of the system can be derived by assembling the motion equations of all elastic shaft elements and discs as follows: (12)

${\kern 0pt} \left\{ \begin{matrix} {M_{de}}{{\ddot u}_x} + \Omega {J_{de}}{{\dot u}_y} + {K_{de}}{u_x} = {Q_{dx}}, \hfill \\ {M_{de}}{{\ddot u}_y} - \Omega {J_{de}}{{\dot u}_x} + {K_{de}}{u_y} = {Q_{dy}}, \hfill \\ \end{matrix} \right.$

where u _x =[x ₁, θ_y ₁, x ₂, θ_y ₂, x ₃, θ_y ₃, x ₄, θ_y ₄]^T, u _y =[y ₁, −θ_x ₁, y ₂, −θ_x ₂, y ₃, −θ_x ₃, y ₄, −θ_x ₄]^T, M _de, J _de, and K _de are the mass matrix, gyroscopic matrix, and stiff matrix, respectively, including all elastic shaft elements and discs. Q _dx and Q _dy are generalized force vectors of discs.

Considering the constrained translational displacements in the x-axis and y-axis of nodes 1 and 4 respectively, the effect of gravity force and seal exciting force, but ignoring tiny rotational displacements, the motion equations of a double disc rotor-seal system can be expressed as follows: (13) $M\ddot u + C\dot u + Ku = Q + G$ , where Q =[Q _d1x+F _sx, Q _d2x+F _sx, Q _d1y+F _sy, Q _d2y+F _sy]^T, Q _d1x=m _d1 e _d1 Ω ²cos(Ωt), Q _d1y=m _d1 e _d1 Ω ²sin(Ωt), Q _d2x=m _d2 e _d2 Ω ²cos(Ωt), Q _d2y=m _d2 e _d2 Ω ²sin(Ωt), u =[x ₂, x ₃, y ₂, y ₃]^T, G =[0, 0, −m _d1 g, −m _d2 g]^T.

Introducing the following non-dimensional transform: $X = x/{c_s},\;Y = y/{c_s},\;t' = \Omega t,$

$\dot x = \Omega {c_s}\dot X,\;\dot y = \Omega {c_s}\dot Y,$

$\ddot x = {\Omega ^2}{c_s}\ddot X,\;\ddot y = {\Omega ^2}{c_s}\ddot Y,$ where c _s is the radial clearance of the seal.

Substituting the expressions into Eqs. (11) and (13), the dimensionless system motion equation can be obtained as (14) $\ddot U + \frac{C}{{\Omega M}}\dot U + \frac{K}{{{\Omega ^2}M}}U = \frac{{Q + G}}{{{c_s}{\Omega ^2}M}}$ , where Q =[m _d1 e _d1 Ω ²cost′+F′_sx, m _d2 e _d2 Ω ²cost′+F′_sx, m _d1 e _d1 Ω ²×sint′+F′_sy, m _d2 e _d2 Ω ²sint′+F′_sy]^T, ${\kern 0pt} \left\{ \begin{matrix} {{F'}_{sx}} \hfill \\ {{F'}_{sy}} \hfill \\ \end{matrix} \right\} = - {c_s}\left[ \begin{matrix} {K_s} - {m_s}{\gamma ^2}{\Omega ^2}\;\;\;\;\;\gamma \Omega {C_s} \hfill \\ \;\;\;\; - \gamma \Omega {C_s}\;\;\;{K_s} - {m_s}{\gamma ^2}{\Omega ^2}\; \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} X \hfill \\ Y \hfill \\ \end{matrix} \right\} - {c_s}\Omega \left[ \begin{matrix} \;\;\;\;{C_s}\;\;\;{\kern 4pt} 2{m_s}\gamma \Omega \hfill \\ - 2{m_s}\gamma \Omega \;\;\;\;\;{C_s}\; \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\dot X} \hfill \\ {\dot Y} \hfill \\ \end{matrix} \right\} - {c_s}{\Omega ^2}\left[ \begin{matrix} {m_s}\;\;\;\;0 \hfill \\ \;0\;\;\;\;{m_s} \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\ddot X} \hfill \\ {\ddot Y} \hfill \\ \end{matrix} \right\}$ , U =[X ₂, X ₃, Y ₂, Y ₃]^T.

Eq. (14) is a second-order differential equation including coefficient matrices, which describes the nonlinear motion of the double disc rotor-seal system. The dimensionless equation of motion is difficult to solve by conventional perturbation methods (Li, ). Therefore, the fourth-order Runge-Kutta method was applied to obtain the numerical response solutions of the system.

4. Numerical results and discussion

Considering the seal parameters K _s, C _s, and γ are nonlinear functions of displacements of the two discs, the motion equation, Eq. (14), shows strongly nonlinear characteristics. The bifurcation diagrams show the variation course of the motion of the system. Time-history diagrams, axis orbit diagrams, Poincaré maps and amplitude spectrums are also presented to analyze the nonlinear vibration of the double disc rotor-seal system. The main parameters and values of the seals, and of the discs and shaft used in the numerical calculation, are shown in Tables 1, 2, and 3, respectively. The experimental coefficients of the seal are listed in Table 4.

Table 1

Main parameters and values of seals

Seal	Radius (m)	Clearance (mm)	Length (m)	Pressure drop (MPa)	Viscosity (Pa∙s)
Left	0.5	0.4	0.06	0.05	0.001
Right	0.3	0.8	0.06	0.05	0.001

Table 2

Main parameters and values of the discs

Disc	Mass (kg)	Dentist (kg/m³)	Mass eccentricity (mm)
Left	120	7830	0.27
Right	180	7830	0.25

Table 3

Main parameters and values of the shaft

Parameter	Value
L ₁ (m)	0.7
L ₂ (m)	0.3
L ₃ (m)	1.0
Diameter (m)	0.5
Dentist (kg/m³)	7830
Elastic modulus (Pa)	2.078×10¹¹

Table 4

Experimental coefficients of seal

Parameter	Value
n ₀	0.079
m ₀	−0.25
γ ₀	0.4
n ₁	2
n ₂	0.5

4.1. Influence of the rotor speed Ω

Rotor speed is an important parameter of rotating machinery which significantly influences the vibration characteristics and dynamic response of rotor-seal systems. To investigate bifurcation, we chose the rotor speed (changing from 500 r/min to 9200 r/min) as abscissa and the dimensionless displacements X ₂ and X ₃ as ordinates. The dimensionless time ranged from 0 to 600π and every period was divided into an average of 256 parts considering the fast Fourier transform (FFT). To ensure accuracy of the results, the preceding half iterative results were discarded. The results showed that the double disc rotor-seal system would exhibit periodic, double-periodic, multi-periodic, and quasi-periodic motion as the rotor speed increases (Fig. 5).

Fig.5
Bifurcation diagrams of left (a) and right (b) discs

The bifurcation diagrams show that the vibration amplitude at a rotor speed of about 1050 r/min is higher because of the occurrence of primary resonance. When the rotor speed is less than 2180 r/min, there is just one solution, which indicates that the system is in period-one motion, that is, the motion of the two discs occurs with the same frequency as the rotor speed. However, once the rotor speed runs over 2180 r/min, the system loses stability and turns into period-two motion. In the range [2320, 2610] r/min, quasi-periodic, period-four, and period-eight motion appear. With increasing rotor speed, the period-two motion appears again until the rotor speed reaches 3120 r/min. Quasi-periodic, multi-periodic (including period-six and period-three motion), and chaotic motion are found in the range [3120, 5290] r/min. Then, the motion of the two discs returns to period-one bifurcation and subsequently changes from period-one motion to quasi-periodic and chaotic motion as the rotor speed increases. Finally, chaotic motion is encountered when the rotor speed is greater than 8310 r/min. The motion status of the two discs is also represented by the time-history diagrams, the axis orbit diagrams, the Poincaré maps and amplitude spectrums (Figs. 6–9, p.45-46).

Fig.6
Numerical analysis of left (a) and right (b) discs at Ω=1200 r/min

Fig.7
Numerical analysis of left (a) and right (b) discs at Ω=4420 r/min

Fig.8
Numerical analysis of left (a) and right (b) discs at Ω=6150 r/min

Fig.9
Numerical analysis of left (a) and right (b) discs at Ω=8650 r/min

Fig. 6 illustrates the dynamic response of the double disc rotor-seal system at Ω=1200 r/min. The system has the same whirl frequency as the frequency of excitation force and has an isolated point in the Poincaré map corresponding to one part of the frequency of the amplitude spectrum, collectively demonstrating that the system exhibits stable synchronous vibration. The axis orbits of the two discs are clearly elliptical, whirling around the points (0.42, −0.77) and (0.19, −0.46), respectively. When the rotor speed increases to Ω=4420 r/min (Fig. 7), the time-history shows the system no longer has stable synchronous vibration. The vibration period becomes three times greater than the period of the excitation force at this rotor speed and the axis orbits begin to diverge. Three isolated points in the Poincaré maps and three isolated frequency parts indicate period-three motion. Fig. 8 shows the numerical analysis results at Ω=6150 r/min, representing the quasi-periodic motion. The axis orbits of the two discs show irregular motion and no obvious period can be found in the time-history diagrams. There is a closed curve in the Poincaré maps and some frequency parts have no common divisor in the amplitude spectrum, clearly indicating that the motion of the system is quasi-periodic. At Ω=8650 r/min, the vibration amplitudes of the two discs are much higher than previously. The axis orbits are becoming completely disorganized. The Poincaré point no longer has finite isolated points or a closed curve, and in the amplitude spectrum diagram, continuous spectrum occurs, collectively demonstrating that the system exhibits chaotic motion.

To study the specific effects of other key parameters on the double disc rotor-seal system, we obtained the bifurcations of both discs, in which the distance of the right disc L ₃, the mass of disc m _d2, the seal clearance c _s2, and the seal clearance L _s2 are labeled as the abscissa, i.e., as variable parameters.

4.2. Influence of the distance of the right disc L ₃

The bifurcation behavior diagrams are presented in Fig. 10 for when the distance of the right disc L ₃ is employed as the variable parameter, ranging from 0.2 m to 1 m. The distance between the left and right discs also varies with L ₃ because the distance of the right disc L ₁ is fixed. The other main parameters are L ₁=0.6 m, m _d2=180 kg, c _s2=0.8 mm, L _s2=0.06 m, and Ω=4000 r/min. Other parameters remain unchanged. The main motion characteristics can be summarized from Fig. 10. As L ₃ increases: multi-periodic motion→quasi-periodic motion→multi-periodic motion→quasi-periodic motion→periodic motion→multi-periodic motion→periodic motion. Note that the system shows stable synchronous vibration when L ₃ ranges from 0.52 m to 0.66 m, which reveals that the stability of a symmetrical double disc rotor-seal system is better than that of an asymmetrical system. When L ₃ is higher than 0.76 m, stable synchronous vibration becomes the main motion, though period-two motion occurs within a narrow range.

Fig.10
Bifurcation diagrams of left (a) and right (b) discs with increasing L ₃

4.3. Influence of the mass of disc m _d2

The bifurcation behavior diagrams are presented in Fig. 11 for when the mass of disc m _d2 is employed as the variable parameter, ranging from 20 kg to 150 kg. The other main parameters are L ₁=0.7 m, L ₃=1 m, c _s2=0.8 mm, L _s2=0.06 m, and Ω=3000 r/min. The main motion characteristics can be summarized as follows. As m _d2 increases: periodic motion→ multi-periodic motion alternating with quasi-periodic motion. Clearly when the mass of disc m _d2 is less than 34 kg, the system shows stable synchronous vibration. But as m _d2 increases, the system not only loses its stability, but the vibration amplitude also becomes larger and the motion state changes frequently in certain zones, ranging from [54.5, 63] to [72.5, 85] kg. The results indicate that a small mass is good for the stability of a double disc rotor-seal system.

Fig.11
Bifurcation diagrams of left (a) and right (b) discs with increasing m _d2

4.4. Influence of the seal clearance c _s2

Fig. 12 shows the bifurcation behavior diagrams for when the seal clearance c _s2 is adopted as the variable parameter, ranging from 0.3 mm to 0.9 mm. The other main parameters are L ₁=0.7 m, L ₃=1 m, m _d2=180 kg, L _s2=0.06 m, and Ω=3000 r/min. The main motion characteristics can be summarized as follows. As c _s2 increases: multi-periodic motion→quasi-periodic motion→multi-periodic motion→periodic motion→multi-periodic motion→quasi-periodic motion→multi-periodic motion. The diagrams indicate that a proper seal clearance c _s2, ranging from 0.376 mm to 0.54 mm, can make the double disc rotor-seal system stable. Quasi-periodic motion also occurs in narrow zones and as c ₂ increases, period-four motion emerges and changes to period-two motion when c _s2=0.617 mm and 0.662 mm, respectively.

Fig.12
Bifurcation diagrams of left (a) and right (b) discs with increasing c _s2

4.5. Influence of the seal length L _s2

Fig. 13 shows the bifurcation behavior diagrams for when the seal length L _s2 is adopted as the variable parameter, ranging from 0.06 m to 0.15 m. The other main parameters are L ₁=0.7 m, L ₃=1 m, m _d2=180 kg, c _s2=0.8 mm, and Ω=3000 r/min. With increasing L _s2, the main motion characteristics can be summarized as follows: period-two motion→quasi-periodic motion→periodic motion. Fig. 13 shows that compared with the previous bifurcation diagrams, the vibration form is single and the system shows period-two motion over a large range as the seal length L _s2 increases, turning to quasi-periodic motion when L _s2 reaches 0.121 m. The oscillating region of the system then starts to decrease and becomes periodic motion when L _s2 is larger than 0.13 m. The results indicate that a larger seal length is beneficial for improvement of the system.

Fig.13
Bifurcation diagrams of left (a) and right (b) discs with increasing L _s2

4.6. Influence of the seal drop pressure ΔP ₂

The bifurcation behavior diagrams are presented in Fig. 14 for when the seal drop pressure ΔP ₂ is employed as the variable parameter, ranging from 0.05 MPa to 0.4 MPa. The other main parameters are L ₁=0.7 m, L ₃=1 m, m _d2=180 kg, L _s2=0.06 m, Ω=3000 r/min, and c _s2=0.8 mm. The main motion characteristics can be summarized as follows. As ΔP ₂ increases: quasi-periodic motion→multi-periodic motion→periodic motion. The diagrams clearly indicate that within the range [0.05, 0.4] MPa, the rotor-seal system becomes more stable with a higher seal drop pressure than with a lower one. Unstable quasi-periodic and multi-periodic motion occurs within the range [0.05, 0.104] MPa and the system becomes stable with increasing ΔP ₂.

Fig.14
Bifurcation diagrams of left (a) and right (b) discs with increasing ΔP ₂

4.7. Comparing linear and nonlinear models

Noah and Sundararajan () discussed the limitations of linearized analysis and the significance of considering nonlinear effects in predicting the dynamic characteristics of rotating systems. They pointed out that nonlinear analysis presents a more realistic representation of a rotating system than linear analysis. Fluid-film forces, including journal bearings, squeeze-film dampers, and annular liquid/gas seals, are highly nonlinear functions of journal displacement and velocity. The linear seal fluid dynamic force model proposed by Childs () was adopted to describe the linear seal force. The linear model can be expressed as follows: (15) ${\kern 0pt} \left\{ \begin{matrix} {F_{sx}} \hfill \\ {F_{sy}} \hfill \\ \end{matrix} \right\} = - \left[ \begin{matrix} \;{K_1}\;\;\;\;{K_2} \hfill \\ - {K_2}\;\;\;{K_1} \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} x \hfill \\ y \hfill \\ \end{matrix} \right\} - \left[ \begin{matrix} \;{C_1}\;\;\;\;\;{C_2} \hfill \\ - {C_2}\;\;\;\;{C_1} \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\dot x} \hfill \\ {\dot y} \hfill \\ \end{matrix} \right\} - \left[ \begin{matrix} {m_s}\;\;\;\;0 \hfill \\ \;0\;\;\;\;{m_s} \hfill \\ \end{matrix} \right]\left\{ \begin{matrix} {\ddot x} \hfill \\ {\ddot y} \hfill \\ \end{matrix} \right\}$ , where K ₁=μ ₃(μ ₀−μ ₂ Ω ² T ²/4), K ₂=μ ₁ μ ₃ ΩT/2, C ₁=μ ₁ μ ₃ T, C ₂=μ ₂ μ ₃ ΩT ², m _s=μ ₁ μ ₃ T ², μ ₀, μ ₁, μ ₂, μ ₃, and T are expressed through the short bearing model (Childs, ).

The bifurcation diagrams of the system using a linear seal force model are shown in Fig. 15. The parameters used in the system are in accord with those of the nonlinear model. The bifurcation diagrams demonstrate that when the rotor speed is less than about 1000 r/min, the vibration amplitudes of the two discs are similar to those of the nonlinear model and the system exhibits stable motion. However, as the rotor speed continues to rise, the vibration amplitudes increase first, and then decrease until the rotor speed reaches 2400 r/min. When the rotor speed runs over 2400 r/min, the system changes from stable period-one motion to unstable motion. The system loses some vibration characteristics compared with the bifurcation diagrams in Fig. 5, which is consistent with the results of Noah and Sundararajan (). Not only are vibration characteristics lost, but also the rub phenomenon emerges. In fact, Child ()’s linear seal model strictly applies only for small motion about the centered position and the nonlinear seal force model can be transformed into a linear seal force model under specific conditions. Thus, Child ()’s linear seal model is an instantiation of the nonlinear seal model of Muszynska (). Although the dynamic coefficients of the linear seal model consider the effects of seal clearance and fluid inertia, they are still far from the actual seal force compared with Muszynska ()’s nonlinear seal model, especially when solving nonlinear problems.

Fig.15
Bifurcation diagrams of left (a) and right (b) discs with a linear model

5. Conclusions

This paper deals with nonlinear vibration characteristics of a double disc rotor-seal system with coupled nonlinear seal forces. The model and the motion equations of the double disc rotor-seal system are established by FEM and the Lagrange equation. Compared with traditional modeling methods, this method is more convenient for studying complex multi-stage rotor-seal systems. The fourth-order Runge-Kutta method is used to solve the motion equations of the system and numerically determine the vibration response of the two discs. The nonlinear dynamic characteristics of the double disc rotor-seal system are represented by bifurcation diagrams, time-history diagrams, axis orbit diagrams, Poincaré maps, and amplitude spectrums, with the rotor speed increasing. The results show that the coupled seal force has great influence on the dynamic stability of the rotor-seal system. The effects of the key parameters, including distance between the two discs, mass of the disc, seal clearance, seal length, and seal drop pressure, are also analyzed. The numerical results show that the system exhibits rich forms of periodic, multi-periodic, quasi-periodic, and chaotic motion, and that a symmetrical disc structure, small disc mass, proper seal clearance, long seal length, and high seal drop pressure are propitious for improving the stability of the system and avoiding severe vibration. Compared with a linear seal model, the nonlinear seal model is more suitable for solving nonlinear vibration problems and analyzing actual situations. These results provide a theoretical foundation for further research on multi-stage rotor-seal systems.

In actual structures of rotating machinery with multi-stage sealing, in addition to the sealing of the impellers, bearings are applied at both ends to support the rotor system. These have a significant impact on the vibration performance of the rotor system. For future research to improve multi-stage rotor-seal systems, further work is needed on the effect of such bearings.

* Project supported by the National Science & Technology Pillar Program during the Twelfth Five-Year Plan Period (No. 2011BAF03B01), China

References

[1] Childs, D.W., 1983. Dynamic analysis of turbulent annular seals based on Hirs’ lubrication equation. Journal of Lubrication Technology, 105(3):429-436.

[2] Childs, D.W., Graviss, M., Rodriguez, L.E., 2007. Influence of groove size on the static and rotordynamic characteristics of short, laminar-flow annular seals. , Joint Tribology Conference of the Society of Tribologists and Lubrication Engineers, American Society Mechanical Engineers, San Antonio, TX. ASME, New York, USA, 422-429. :422-429.

[3] Ding, Q., Cooper, J.E., Leung, A.Y.T., 2002. Hopf bifurcation analysis of a rotor/seal system. Journal of Sound and Vibration, 252(5):817-833.

[4] Fei, Z.X., Tong, S.G., Wei, C., 2013. Investigation of the dynamic characteristics of a dual rotor system and its start-up simulation based on finite element method. Journal of Zhejiang University-SCIENCE A (Applied Physics & Engineering), 14(4):268-280.

[5] Hua, J., Swaddiwudhipong, S., Liu, Z.S., 2005. Numerical analysis of nonlinear rotor-seal system. Journal of Sound and Vibration, 283(3-5):525-542.

[6] Kaneko, S., Ikeda, T., Saito, T., 2003. Experimental study on static and dynamic characteristics of liquid annular convergent-tapered damper seals with honeycomb roughness pattern. Journal of Tribology, 125(3):592-599.

[7] Li, S.T., Xu, Q.Y., Zhang, X.L., 2007. Nonlinear dynamic behaviors of a rotor-labyrinth seal system. Nonlinear Dynamics, 47(7):321-329.

[8] Li, W., Yang, Y., Sheng, D.R., 2011. Nonlinear dynamic analysis of a rotor/bearing/seal system. Journal of Zhejiang University-SCIENCE A (Applied Physics & Engineering), 12(1):46-55.

[9] Megerle, B., Rice, T.S., McBean, I., 2013. Numerical and experimental investigation of the aerodynamic excitation of a model low-pressure steam turbine stage operating under low volume flow. Journal of Engineering for Gas Turbines and Power, 135(1):012602

[10] Muszynska, A., 1988. Improvements in lightly loaded rotor/bearing and rotor/seal models. Journal of Vibration, Acoustics, Stress, and Reliability in Design, 110(2):129-136.

[11] Muszynska, A., Bently, D.E., 1990. Frequency-swept rotating input perturbation techniques and identification of the fluid force models in rotor/bearing/seal systems and fluid handling machines. Journal of Sound and Vibration, 143(1):103-124.

[12] Nelson, H.D., McVaugh, J.M., 1976. The dynamics of rotor-bearing systems using finite elements. Journal of Engineering for Industry, 98(2):593-600.

[13] Noah, S.T., Sundararajan, P., 1995. Significance of considering nonlinear effects in predicting the dynamic behavior of rotating machinery. Journal of Vibration and Control, 1(4):431-458.

[14] Smalley, A.J., Camatti, M., Childs, D.W., 2006. Dynamic characteristics of the diverging taper honeycomb-stator seal. Journal of Turbomachinery, 128(4):717-724.

[15] Wang, W.Z., Liu, Y.Z., Meng, G., 2009. A nonlinear model of flow-structure interaction between steam leakage through labyrinth seal and the whirling rotor. Journal of Mechanical Science and Technology, 23(12):3302-3315.

[16] Wang, W.Z., Liu, Y.Z., Meng, G., 2009. Nonlinear analysis of orbital motion of a rotor subject to leakage air flow through an interlocking seal. Journal of Fluids and Structures, 25(5):751-765.

[17] Wang, Y.F., Wang, X.Y., 2010. Nonlinear vibration analysis for a Jeffcott rotor with seal and air-film bearing excitations. Mathematical Problems in Engineering, 2010:657361

[18] Yuan, Z., Chu, F., Hao, R., 2007. Clearance-excitation force of shrouded turbine rotor accounting for pitching motion. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 221(2):187-194.

[19] Zhang, W., 1990. The Theoretical Base of Rotordynamic. (in Chinese), Science Press,Beijing :

[20] Zhong, Y.E., 1987. Rotordynamics. (in Chinese), Tsinghua University Press,Beijing :

Open peer comments: Debate/Discuss/Question/Opinion

<1>

Similar articles