基于有限元法的转子轴承系统非线性特性研究摘要针对典型的转子轴承系统构造了一个复杂多因素并且能够比较真实地反映实际系统的非线性系统模型。

采用有限元方法将其离散化分为圆盘、轴段和轴承座等单元，并对各单元作了详细的动力分析，当考虑油膜力耦合作用时，广义力的求解引用了瑞利耗散函数，推出了油膜粘性阻尼力的非线性因素，再由拉格朗日方程得出系统的运动微分方程。

最后关键词：陀螺力矩油膜力转子轴承系统有限元Finite element method based on nonlinear characteristics of rotor bearing Abstract A typical rotor-bearing system for a complex multi-factor structure and the ability to truly reflect the actual system of nonlinear system model. Finite element method to the disc is divided into discrete, such as shafts and bearing units, each unit made a detailed and dynamic analysis, when considering the coupling of oil film force, the generalized Rayleigh power dissipation of the solution quoted function, introduced the film's nonlinear viscous damping factor, then the Lagrange equations derived differential equations of motion. Finally,Key words: oil film force gyroscopic element rotor-bearing system1 引言转子系统在机械、动力、航空航天等领域有着广泛的应用，是机器设备的重要组成部分，随着旋转机械向高速、大功率的方向发展，在旋转机械中常常会出现非线性动力学现象(例如:跳跃、分岔和混沌等)，其对设备的运行构成了严重的威胁。

因此转子动力系统的稳定性成为人们日益关注的问题。

轴承一转子系统是一个复杂的非线性动力系统。

文献[1]研究了非线性轴承-转子系统运用时间有限元法对一个径向游隙的轴承模型与挠性轴的有限元模型求解出了系统的不平衡响应。

文献[2]就600MW汽轮机组转子-轴承系统，建立了系统的运动方程和转子模型，采用有限元分析软件ANSYS 进行模态分析，计算汽轮机转子轴承系统的固有频率和临界转速，分析了转子的特性。

文献[7,8]研究了转子动力学中轴系弯扭耦合的一些非线性动力特性。

本文采用有限元法将转子轴承系统划分了3大单元，综合考虑了系统中存在的油膜力、陀螺力、不平衡力等严重的非线性激励源，建立了比较复杂的数学模型。

最后采用数值分析法求解系统的运动微分方程，并给出了仿真实验。

2 转子轴承系统动力学模型一个典型的转子-轴承系统通常可以沿轴线把转子系统划分为圆盘、轴段和轴承座等单元[3]。

各单元间彼此在结点处连结。

这些结点通常是选在圆盘中心，轴颈中心以及轴线的某些位置上，并按顺序编号（如图1）。

图1转子轴承系统以轴承座中心线为s轴，建立固定坐标系oxys。

转子轴的任一横截面位置可由如下两个位移向量表示，其中x、y为轴心坐标，x yθθ、为截面的偏转角，以及自转角ϕ表示。

2.1 圆盘设圆盘轴心与重心重合，圆盘的广义坐标是其轴心结点的位移向量，{}1,Td yu xθ⎡⎤=⎣⎦和{}[]2,Td xu yθ=-。

oξηζ'以轴心结点为原点，固结在圆盘上的动坐标系（如图2）1图2 圆盘上的动坐标系引入广义坐标并略去高阶小量后得圆盘动能：{}[]{}{}[]{}{}[]{}1122212112212=++Ω+ΩTTd d d d d Td d d d p T uM u u M uu J u J （1）式中 m 为圆盘质量，J d 为赤道转动惯量，J p 为极转动惯量， []00d d mM J ⎡⎤=⎢⎥⎣⎦[]000p J J ⎡⎤=⎢⎥⎣⎦由质量不平衡所对应的广义力近似为{}{}2122Q cos sin 00Q cos sin 00u du de e m t t e e m t t ξηηξ⎫⎛-⎫⎧⎫⎧⎫=ΩΩ+Ω⎪⎨⎬⎨⎬ ⎪⎩⎭⎩⎭⎝⎭⎪⎬⎛⎫⎧⎫⎧⎫⎪=ΩΩ+Ω⎨⎬⎨⎬ ⎪⎪⎩⎭⎩⎭⎝⎭⎭（2）式中e ξη、e 为质量偏心距在坐标系o ξηζ'下的坐标.拉格朗日方程如下式所示：ui ii i F q Vq T q T dt d =∂∂+∂∂-⎪⎪⎭⎫ ⎝⎛'∂∂（3）其中F ui 为保守力。

由拉格朗日方程可得刚性圆盘的运动微分方程[][]1112220000u d d d d u d d d d M u uJ Q M u uJ Q ⎡⎤Ω⎧⎫⎡⎤⎧⎫⎧⎫+=⎨⎬⎨⎬⎨⎬⎢⎥⎢⎥-Ω⎣⎦⎩⎭⎩⎭⎩⎭⎣⎦ （4）其中[]d M 为圆盘的惯性力矩，[][]g M J =Ω为陀螺力矩，Ω为转子自转角速度2.2 轴承轴承座简化成图3示单元，轴承座中心坐标是b b x y 、，轴颈中心与其重合。

图3轴承座单元若不计阻尼影响且认为支承是各向同性的，则有bxy byx yx xy k k k k ====0bxy byx yx xy c c c c ====xx yyk k =bxx byy k k =bxx byy c c =，由此推得轴承座的运动方程是：0000000b b bxbxxb b by byy b b ux bxx xxb b ubyy yy b b y M c x x M c y y Q k k x x k k y y Q ⎡⎤⎡⎤⎧⎫⎧⎫+⎨⎬⎨⎬⎢⎥⎢⎥⎩⎭⎩⎭⎣⎦⎣⎦⎧⎫⎡⎤⎡⎤⎧⎫⎧⎫⎪⎪+==⎨⎬⎨⎬⎨⎬⎢⎥⎢⎥⎩⎭⎩⎭⎣⎦⎣⎦⎪⎪⎩⎭（5） 2.3 轴单元分析采用铁摩辛柯梁模型，计入轴的剪切变形，分别用形状函数表示轴单元的动能和势能。

2.3.1形状函数如图4所示，()11,,,j j j j y y ψψ--为轴单元在其左右两端面上的位移、转角。

()j y ψ⎛⎫ ⎪ ⎪11j j y ψ--⎛⎫⎪ ⎪,,,E I A ρ图4 第j 个轴单元图5弹性轴段单元该单元上任意一点s 处的位移或挠度表征为端面坐标的函数，()()111234j j j j y y s N N N N y ψψ--⎛⎫ ⎪⎪= ⎪ ⎪ ⎪⎝⎭（6）记为 ()j y s NY = 其中()1,2,3,4i N i =称为形状函数。

2.3.2 轴的动能如图5所示为一弹性轴段单元，该单元的广义坐标是两端结点的位移，即{}{}[]12,,,,,,Ts A yA B yB Ts A xA B xB u x x u y y θθθθ⎫⎡⎤=⎪⎣⎦⎬⎪=--⎭（7）单元内任一截面的位移,,,y x x y θθ是位置s 和时间t 的函数。

单元的结点位移可用形状函数来表示：()[]{}1,s x s t N u =()[]{}2,s y s t N u =又[]{}1ys x N u sθ∂'==∂[]{}2xs y N u sθ∂'-==∂得轴单元的动能11221122/21,/22,2c s s T Ts s s c s s s s T Ts s ss M J uT u u M J uJ u u u J u +⎡⎤⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦+⎣⎦⎣⎦⎡⎤⎡⎤Ω⎡⎤+⎢⎥⎢⎥⎣⎦-⎣⎦⎣⎦（8） 式中[][][]123434621250(1)Tc m m m m m m m lAl M A N N ds m m m ρρφ⎛⎫⎪-⎪==⎪-+ ⎪⎝⎭⎰对称[][][]()78789810227892021Ts lJ I N N dsm m m m m m m A lI m m A l m ρρφ''=-⎛⎫⎪- ⎪=⎪-+ ⎪⎝⎭⎰对称2113735103m φφ=++22111121012024m l φφ⎛⎫=++ ⎪⎝⎭239370106m φφ=++241334204024m φφ⎛⎫=-++ ⎪⎝⎭225110560120m l φφ⎛⎫=++ ⎪⎝⎭226114060120m l φφ⎛⎫=-++ ⎪⎝⎭ 765m =81102m l φ⎛⎫=- ⎪⎝⎭22921563m l φφ⎛⎫=++ ⎪⎝⎭221013066m l φφ⎛⎫=--+ ⎪⎝⎭212EIAG lφμ=2.3.3 轴的势能势能只是位置的函数，她代表着内力和外力所做的总功，轴单元的势能也可用形状函数和其端面坐标表示，由功能互等原理可得： {}[]{}{}[]{}11221122TTs s s s s s s V u K u u K u =+（9） 式中[][][]2232126126(4)62(2)126(1)(4)Ts l K EI N N dslll l l EI l l l φφφφ''''=-⎛⎫ ⎪+-- ⎪=⎪-+ ⎪+⎝⎭⎰对称由Lagrange 方程可推出轴段单元的运动方程：1122111222/2/2c s s s s c s s ss us s u s s M J u J uM J u J uu K Q u K Q +⎡⎤⎡⎤⎡⎤⎡⎤+Ω⎢⎥⎢⎥⎢⎥⎢⎥+-⎣⎦⎣⎦⎣⎦⎣⎦⎧⎫⎡⎤⎡⎤+=⎨⎬⎢⎥⎢⎥⎣⎦⎣⎦⎩⎭（10）其中{}[]21cos sin 0Tus lQ N e t e t dsξημ⎡⎤=ΩΩ-Ω⎣⎦⎰ {}[]22cos sin 0TuslQ N e t e t dsηξμ⎡⎤=ΩΩ+Ω⎣⎦⎰3 转子轴承系统的运动方程对于具有N 个结点，N1个轴承支承和N2个圆盘，其间用N-1个轴段连接而成的转子系统，则系统的位移向量是{}11112222,,[,,,,,,,,...,,]y x y Tx N yN N xN U x y x y x y θθθθθθ=---（11）综合各圆盘、轴段单元及轴承支承的运动方程（4）式、(10)式及（5）式，可得转子系统的运动方程：[]{}[]{}[]{}{}M UC U K U Q ++= （12）此时[M]、[C]、[K]都是半带宽为8的稀疏带状矩阵，故求解比较方便。

4 计算仿真及分析用标准四阶龙格库塔算法对式(12)进行数值积分,可求得系统的自由振动响应，在分析时,计算结果以分岔图、poincare 映射图和频谱图等形式给出。