不同坐标系下滑动轴承动特性系数计算方法的相互转换*

王蕴馨 马金奎 陈淑江 路长厚 李 佳 刘志颖

(山东大学机械工程学院,高效洁净机械制造教育部重点实验室 山东济南 250061)

目前滑动轴承动力特性研究常用的计算方法有差分法和偏导数法。为避免两次差分近似计算互相叠加的计算误差[3]，本文作者以圆轴承为对象，采用有限差分法求解雷诺方程，在计算油膜压力过程中，讨论了雷诺边界条件的实现方法；同时根据偏导数法计算扰动压力并给出了其分布，从而得到滑动轴承的动特性系数。

在研究滑动轴承流体润滑时，油膜厚度和动特性系数的表示通常有2种坐标系，即极坐标系[4]和直角坐标系[5]，本文作者分别建立了2种坐标系下的数值分析模型以及二者之间的转换矩阵，经验证可获得相近的结果。文中的研究完善了经典轴承理论的坐标转换方法，对滑动轴承的研究有一定的指导作用。

1 油膜压力分布

1.1 油膜厚度与雷诺方程的量纲一化形式

图1为油膜厚度计算简图，Ob为轴承中心，Oj为轴的中心；点A为沿y轴正方向逆时针转动θ后的计算点;φ为偏位角，φ=θ-φ；e为偏心距；ω为轴的转速。设c为半径间隙，可得计算点A的膜厚表达式为

(1)

图1 油膜厚度计算简图Fig 1 A sketch map of oil film thickness calculation

根据滑动轴承的流体润滑理论，将雷诺方程量纲一化，得

(2)

1.2 雷诺方程的求解及边界条件

采用有限差分法求解雷诺方程[6-7]，将油膜划分为m×n个网格，沿φ方向划分m等分并用i编号，沿λ方向划分n等分并用j编号，节点位置用(i,j)二维编号表示。利用差商取代雷诺方程中的一阶偏导数，整理可得

(3)

其中：

Di,j=Ai,j+Bi,j+Ci,j+Ci,j

Ei,j=3(Hi+1,j-Hi,j)Δφ

计算油膜压力时，边界条件的合理采用是影响滑动轴承润滑性能计算精度的重要因素。文中引入雷诺边界条件[8-10]：在油膜破裂边界上，油膜的压力和压力的一阶导数均为0。即：

P=0，∂P/∂φ=0，φ∈[0,2π]

如图2所示，实现雷诺边界条件的具体方法是：计算油膜压力时，由起始点向终止点方向逐点迭代。先假设内部各点压力为0，运用有限差分法依次计算各点的压力值，如果算出某点压力值为负数，则判定此点出现油膜破裂并取为0，这样就得到了油膜压力的第1次近似值。利用该近似值再次从第一个点重新计算压力值，如此反复迭代直至满足收敛误差。可以看出，油膜开始破裂的位置随着迭代次数的增加逐渐向下游推移，则整个压力曲线随之增高，于是破裂边界和压力分布会逐渐逼近实际破裂边界和雷诺边界条件的压力分布。

图2 渐变的压力曲线下游Fig 2 Downstream of gradual pressure curves

从文献[11-12]中取一组轴承参数作示例，如表1所示。所得量纲一油膜压力分布如图3所示。

表1 轴承参数Table 1 Parameters of a bearing

图3 量纲一油膜压力分布Fig 3 Distribution of dimensionless oil film pressure

可以看出，滑动轴承的量纲一油膜压力的三维分布近似抛物面分布。在0≤φ≤π的区域，量纲一油膜压力在某一段逐渐增大到最大压力值，之后急剧下降，在φ>π的某一区域降为0。这与经典滑动轴承润滑理论一致[13-14]。

2 偏导数法求解滑动轴承动力特性系数

油膜具有非线性力学特性，其刚度和阻尼会影响系统的运动稳定性[15]。油膜刚度和阻尼系数的2种坐标系下的求解模型如图4所示。

图4 滑动轴承动特性系数计算模型Fig 4 A calculation model of dynamic characteristic coefficients of sliding bearing

图4中，Ob为轴承中心，O0为静平衡位置,Oj为轴的瞬时中心，极坐标系为εObφ;直角坐标系为xOby。轴心做微小运动时的量纲一油膜力为

(4)

(5)

2.1 极坐标系求解模型

(6)

定义油膜刚度系数为单位位移所引起的油膜力增量，即

定义油膜阻尼系数为单位速度所引起的油膜力增量，即

其中，下标0表示在静平衡位置处求导；各系数的第1个下标代表力的方向，第2个下标代表位移或速度的方向(下同)。

(7)

将扰动压力在完整油膜区内积分，可得到极坐标系中的量纲一刚度系数和阻尼系数，即

(8)

(9)

计算时，扰动压力的边界条件是：在完整油膜区的全部周边上均等于0。根据表1数据及上述分析，可得到滑动轴承量纲一扰动压力周向分布，如图5所示。

图5 量纲一扰动压力周向分布Fig 5 Circumferential distribution of dimensionless disturbed pressure (a)circumferential distribution of Pε;(b)circumferential distribution of Pφ;(c)circumferential distribution of

2.2 直角坐标系求解模型

同理，将式(5)在静平衡位置展开为Taylor级数，仅保留一阶微量，可得到

(10)

同样地，可以定义

(11)

扰动压力的边界条件不变，积分后可得到直角坐标系中的量纲一刚度系数和阻尼系数，即

(12)

(13)

3 坐标转换矩阵

3.1 坐标转换矩阵推导

(14)

(15)

(16)

已知油膜力可表示为

因此

根据系数相等可得2种数值分析模型的转换方法为

(17)

或

(18)

经典轴承理论[1]中的转换矩阵如式(19)所示，由于结论有所差异，接下来给出验证分析。

(19)

3.2 验证分析

首先确定边界条件，求解压力分布，然后根据偏导数法得到的扰动压力雷诺方程的量纲一化形式计算各扰动压力，再经积分求出2种计算模型的量纲一动特性系数，最后通过坐标转换矩阵换算后的结果，如表2所示。

表2 量纲一动特性系数(l/d=0.5,φ=30°)Table 2 Dimensionless dynamic characteristic coefficients(l/d=0.5,φ=30°)

2种计算模型得到的量纲一动特性系数近似相等，最大误差为6.05%，其中，交叉阻尼系数近似相等，这有力地论证了上述坐标转换方法的正确性。

4 结论

(1)分别建立极坐标系和直角坐标系下滑动轴承动特性系数的计算模型，2种模型均可求解出滑动轴承油膜压力和扰动压力分布及其动特性系数，同时给出了二者之间的坐标转换矩阵，并进行了验证，说明了方法的正确性。

(2)通过对压力曲线下游边界变化情况的分析，结合迭代算法，研究了雷诺边界条件的实现方法：先假设内部各点压力为0，运用有限差分法和负压置零的条件逐点计算出各点的第1次压力近似值，利用该近似值再次从第一个点重新计算，如此反复迭代直至满足收敛误差，从而得到满足条件的破裂边界和压力分布。