滑动轴承-转子系统Riccati-Newmark加速度传递矩阵法

第一作者毛文贵女，博士生，副教授，1975年1月生

通信作者韩旭男，博士,教授，博士生导师，1968年6月生

滑动轴承-转子系统Riccati-Newmark加速度传递矩阵法

毛文贵1,2，韩旭1，刘桂萍1

(1. 湖南大学机械与运载工程学院,长沙410082；2. 湖南工程学院机械工程学院,湖南湘潭411101)

摘要：为克服传递矩阵法数值不稳定及非线性瞬态响应分析中滑动轴承矩阵建立困难问题，提出Riccati-Newmark加速度传递矩阵法。借助Newmark加速度法建立传递矩阵，采用Taylor级数预估滑动轴承轴心下一时刻位移、速度建立滑动轴承矩阵；据边界条件用Riccati传递矩阵法求滑动轴承-转子系统非线性瞬态响应，提高数值稳定性。以单圆盘转子系统为例，与传统轴颈下一时刻位移、速度近似线性扰动处理的瞬态响应对比分析，验证此方法的有效性；讨论不同转速下线性、非线性油膜力的瞬态轨迹。

关键词：非线性转子系统；瞬态响应；Newmark加速度法；Riccati传递矩阵法

基金项目：国家自然科学基金11202073;高等学校博士学科点专项科研基金(20130161130001);“高档数控机床与基础制造装备”科技重大专项2012(2012ZX04002-091)

收稿日期：2014-08-07修改稿收到日期：2014-09-25

中图分类号:TH113 .3,1;O322

文献标志码：A

DOI:10.13465/j.cnki.jvs.2015.20.014

Abstract:In order to eliminate the numerical instability of transfer matrix method and to build a transfer matrix of nonlinear elements( bearings), the Riccati transfer matrix method was extended to the transient analysis of nonlinear rotor-bearing systems. In the method, the transfer matrix was obtained with the aid of the Newmark acceleration formulation, the deflections and velocities at the stations containing nonlinear element (bearings) were predicted by Taylor series, and then the deflections, velocities and accelerations at all stations were solved by using the Riccati transfer matrix method combined with the Newmark acceleration formulation integration according to the boundary conditions. An example of single disc rotor system was presented and the results were compared with those by a transient analysis considering linear perturbation and linear oil film force to verify the effectiveness of the proposed method.

Riccati transfer matrix method combined with Newmark acceleration formulation integration for analysing sliding bearings and rotor system

MAOWen-gui1,2,HANXu1,LIUGui-ping1(1. HNU College of Mechanical and Vehicle Engineering, Hunan University, Changsha 410082, China；2. College of Mechanical Engineering, Hunan Institute of Engineering, Xiangtan 411101, China)

Key words:nonlinear rotor system; transient response; Newmark acceleration method; Riccati transfer matrix method

滑动轴承作为弹性支承可降低轴系结构刚度，使轴系在较宽速度、载荷范围内无磨损工作，具有高回转精度。轴系响应作为研究转子动力学的主要内容，其准确性较大程度上取决于油膜力处理。由于油膜力的强非线性，滑动轴承支承轴系为非线性系统。传递矩阵法[1-2]具有矩阵阶数不随轴系自由度增加而增加、编程简单等优点，最适用转子类链式系统，多用其分析研究线性转子系统动力响应问题，而对非线性转子系统动力响应分析研究有限。对非线性转子系统，滑动轴承传递矩阵只有获得非线性油膜力方能真正建立，非线性油膜力与轴心位移、速度有关，而非线性系统动力响应未知，因此获得滑动轴承传递矩阵较困难。加之滑动轴承强非线性精确表达式不易获得，从而阻碍滑动轴承-转子系统非线性研究进展。转子系统瞬态响应分析中关于轴承油膜力处理有两种方式：线性扰动即对下一时刻轴颈位移用上一时刻位移代替；将线性化油膜力用4个刚度、阻尼系数线性化表示。李苹等[3]采用Riccati传递矩阵与Newmark数值积分法研究滑动轴承-转子系统非线性瞬态响应。考虑油膜非线性，通过CFD软件模拟[4-5]，但对轴心位移、速度用线性扰动近似。线性油膜力、线性扰动均不能真实反映轴系强非线性。用传递矩阵法研究轴系非线性瞬态响应时传递矩阵中含材料、结构及时间等参数，由于一般轴系中材料弹性模量加大而瞬态响应分析中时间间隔较小，采用力、位移状态变量表示传递矩阵易出现较0(1)量级大得多、小得多的数，会引起新病态矩阵[6]。

对此，本文提出非线性扰动(Taylor级数预估)及非线性油膜力的Riccati-Newmark加速度传递矩阵法。用力、加速度状态变量避免数值不稳定性[7]。为获得非线性系统准确的油膜响应，通过Taylor级数预估滑动轴承非线性元件下一时刻位移、速度，借助Newmark加速度法建立传递矩阵，用Riccati传递矩阵法求滑动轴承-转子系统非线性瞬态响应，获得滑动轴承下一时刻位移及速度,通过精度控制误差获得滑动轴承位移、速度。以单圆盘弹性转子系统为例，分别采用线性油膜力、非线性油膜力，线性扰动、非线性扰动对转子系统瞬态响应比较分析，验证本文方法的有效性。

1Newmark加速度法

求解非线性运动微分方程通常用数值解法，如Runge-Kutta法、Newmark法[8]等。前者收敛速度慢，易发散且积分步长对结果精度及正确性影响较大。后者当参数选择合适时积分无条件稳定，结果与步长大小无关。本文用后者求解非线性瞬态响应。

对已知上一步位移、速度、加速度及求解下一步位移、速度、加速度Newmark法进行变换，获得Newmark加速度法，即

(1)

2改进的轴系各元件传递矩阵

传递矩阵法分析轴系结构时将轴系离散为无质量轴段与集总圆盘(F=0)、有支承集总圆盘(F≠0)等单元，见图1。模块中m，Jp，JD为集总圆盘节点处质量、极转动惯量、直径转动惯量；D，d，L为无质量轴段外径、内径及跨度。

图1　转子轴承系统模块化图 Fig.1 bearings-rotor systems modular figure

每种单元均有传递矩阵。Riccati传递矩阵法为提高传递矩阵数值稳定性，将状态向量中元素分为两部分，即一部分由对应初始截面状态向量中具有零值元素组成，另部份由其余互补元素组成。状态向量为

q=(My,Qy,Mz,Qz⋮θy，Y,θz,Z)T=(f⋮e)T

式中：Mz，My为弯矩；Qz，Qy为剪力；Z，Y为位移；θz，θy为转角。

用力、位移状态变量传递矩阵法研究轴系非线性瞬态响应时矩阵中含材料、结构及时间参数，由于弹性模量较大、时间间隔较小，矩阵中出现较0(1)量级大得多、小得多的数会引起新病态矩阵。本文用式(1)改造力、位移状态变量表示的传递矩阵，获得用力、加速度状态变量表示的传递矩阵为

2.1集总圆盘点传递矩阵

由达郞伯尔原理可导出集中圆盘运动微分方程为

(2)

式中：ω为自转角速度；CD为集总圆盘处外阻尼系数；g为重力加速度；eμ为偏心矩。

用式(1)改造式(2)，得用力、加速度状态变量表示的集总圆盘点传递矩阵为

(3)

式中：D1，D2为集总圆盘点传递矩阵，即

式中：Δt为时间间隔。

2.2无质量轴段场传递矩阵

由文献[9]可知无质量轴段y方向左右两端状态变量关系为

(4)

动力学方程为

(5)

式中：K为轴段弹性刚度。

无质量轴段z方向左右两端状态变量关系与动力学方程类似，仅g=零。用式(1)、(5)改造式(4)，获得用力、加速度状态变量表示的无质量轴段场传递矩阵为

(6)

式中：B1，B2为无质量轴段场传递矩阵，即

式中：

V=6EJ/aGAL2为考虑轴段的剪切影响，a为截面系数(实心轴0.886，空心轴2/3)，G为剪切弹性模量，A为截面积

2.3滑动轴承点传递矩阵

滑动轴承类似有支承的集总圆盘(F≠0)，动力学方程为

(7)

由于式(7)中含非线性项，因此将滑动轴承位移、速度及非线性项作为已知量，加速度作为传递变量，得

(8)

式中：C1，C2为滑动轴承点传递矩阵，即

3轴系响应Riccati法

对每一元件(集总圆盘、无质量轴段及滑动轴承)由式(3)、(6)、(8)获得相应传递矩阵，集总圆盘及滑动轴承作为节点处理，与其右边轴段组成典型单元。该单元传递矩阵为

(9)

引入Riccati变换[10]，建立力状态变量及加速度状态变量之关系为

(10)

由式(9)、(10)变换得

(11)

(12)

(13)

(P)j+1=(U11P+Ff)j-(S)j+1(U21P+Fe)j

(14)

(Q)j=-(T)j(U21P+F1e)j。

4Riccati-Newmark加速度传递矩阵法

对线性轴系可据边界及初始位移、速度条件求出各节点初始加速度，再通过轴系响应求解方法计算下一时刻各节点加速度，由式(1)获得各节点位移及速度。而对非线性轴系含非线性元件节点的传递矩(如滑动轴承，式(8))，只在非线性油膜力已知后其节点传递矩阵才真正建立。因此用轴系响应Riccati法计算加速度前须已知非线性项，否则无法计算。对此，本文采用Taylor级数预估非线性元件的下一时刻t+Δt的位移、速度，计算式为

(15)

本文所提Riccati-Newmark加速度传递矩阵法流程见图2。流程中按式(16)修正t+Δt时刻位移及速度，即

(16)

式中：ε为迭代收敛需满足的精度。

采用图2流程获得t+Δt时刻各节点位移、速度。重复可得t+2Δt，t+3Δt，….等各时刻各节点位移、速度及加速度响应。

图2　Riccati-Newmark加速度传递矩阵法 Fig.2 Riccati-Newmark transfer matri acceleration formulation integration method

5数值例子

采用单圆盘转子系统模型：轴承间隙h0=100 um，轴承宽Lc=28.7mm，转轴跨度L=2 000 mm，直径d=57.4 mm；集总圆盘直径D=100 mm, 质量为5 kg；圆盘处偏心距eu=1 mm，轴段质量密度7.80 g/cm3，弹性模量E=211 GPa，泊松比0.3。油膜动力粘度0.016 8 Pa.s。轴段划分72段，73个节点，两轴承分别作用于第6、第66节点，圆盘作用于第36节点；滑动轴承用capone圆轴承模型。采用本文提出的方法对滑动轴承-转子系统进行非线性瞬态响应分析。

在低速1 200 r/min下用本文方法及Runge-Kutta法(ode45位移法)获得轴承瞬态、稳态轨迹见图3。由图3看出，两方法瞬态解完全吻合，稳态轨迹基本吻合。由此验证本文方法的正解性。Runge-Kutta法为对下一时刻轴承处轴颈位移、速度采用线性扰动方式，线性扰动最后稳态下位移较本文方法偏大，用线性扰动处理轴心下一时刻位移相对安全。

针对Taylor级数预估轴心下一时刻位移的非线扰动方法，本文分别对比分析不同转速下线性、非线性油膜力作用的轴心轨迹，结果见图4。由图4看出，较低转速下，开始线性油膜力Y向(重力作用方向)响应较大，稳定后响应范围也较大。因此小扰动下工程中一般用8个系数处理非线性油膜力；但高转速时非线性油膜力响应幅值较线性大。

分析图3、图4结果可知，小扰动下采用线性扰动响应较非线性扰动响应大。用线性扰动模拟下一时刻轴颈处位移、速度瞬态分析设计偏于安全；但线性扰动化处理对高精密结构设计较粗糙。高速下油膜力非线性特性对轴系瞬态响应有一定影响，采用非线性油膜力及非线性扰动的瞬态响应幅值更大。即转速达一定范围时不能视油膜力为线性。

图3 两种方法轴承轴心轨迹Fig.3Orbitofshaftcenterbytwomethods图4 轴承轴心轨迹Fig.4Orbitofshaftcenter

6结论

(1)本文所提非线性滑动轴承-转子系统瞬态响应新方法综合Riccati传递矩阵法与Newmark加速法优点，改进传递矩阵的稳定性。通过与传统瞬态响应分析的轴心下一时刻位移近似扰动所得结果对比，验证该方法的准确性。

(2)通过引入Taylor级数预估，解决非线性元件滑动轴承传递矩阵建立；通过对比不同转速非线性扰动下线性油膜力、非线性油膜力产生的瞬态响应知，高速下油膜力线性处理不能诠释轴系非线性现象。

参考文献

[1]何成兵,顾煜炯.增量传递矩阵法及其在轴系弯扭耦合振动中的应用[J]. 振动工程学报,2006,19(2):119-226.

HE Cheng-bing,GU Yu-jiong.Increment transfer matrix method and its application in the coupled flexural and torsional vibrations analysis[J]. Journal of Vibration Engineering, 2006, 19(2):119-226.

[2]高浩鹏,黄映云,刘鹏.引入传递矩阵法的复杂多体系统连接件建模方法研究[J].振动与冲击,2012,31(6):52-55.

GAO Hao-peng,HUANG Ying-yun,LIU Peng. Modeling of linkers in complex multibody system with transfer matrix method[J].Journal of Vibration and Shock,2012,31(6):52-55.

[3]李苹,窦海波,王正.水轮发电机组主轴系统的建模及非线性瞬态响应[J].清华大学学报:自然科学版,1998,38(6)：123-128.

LI Ping,DOU Hai-bo,WANG Zheng. Modeling for rotor-bearing system of hydroelectric machines and its non-linear transient response[J].Journal of Tsinghua University: Sci. &Tech.,1998,38(6):123-128.

[4]Liu H P, Xu H, Ellison P J,et al.Application of computa-tional fluid dynamica and fluid-structure interaction method to the lubrication study of a rotor-bearing system[J].Tribol Letters,2012,38(3):325-336.

[5]Li qiang,Yu gui-chang,Liu shu-lian. Application of computational fluid dynamics and fluid structure interaction techniques for calculating the 3D transient flow of journal bearings coupled with rotor systems[J].Chinese Journal of Mechanical Engineering,2013,25(5):926-931.

[6]顾致平,刘永寿.非线性转子系统中的传递矩阵技术[M].北京:科学出版社, 2010.

[7]王正. Riccati传递矩阵法的奇点及其消除方法[J]. 振动与冲击, 1987,6(2): 74-78.

WANG Zheng. Singular point and elimination method of riccati transfer matrix method[J]. Journal of Vibration and Shock, 1987,6 (2):74-78.

[8]徐斌. Matlab有限元结构动力学分析与工程应用[M].北京:清华大学出版社, 2010.

[9]闻邦椿,顾家柳,夏松波,等. 高等转子动力学[M]. 北京：机械工业出版社, 2000.

[10]毛文贵,刘桂萍,韩旭,等.一种高效的电主轴系统复频率计算方法[J].振动与冲击,2014,33(19):164-168.

MAO Wen-gui,LIU Gui-ping,HAN Xu,et al.An efficent method for calculating complex frequency of a motorized spindle[J].Journal of Vibration and Shock, 2014, 33(19): 164-168.