动压径向滑动轴承油膜力场的模拟与数值计算研究

张一磊，徐武彬，李 冰，倪 壮

（广西科技大学机械工程学院，广西 柳州 545006）

1 引言

动压径向滑动轴承由于其结构简单、制造方便、能承受重载及冲击等优点被广泛应用在大型燃气轮机、机床、发电机、电动机等旋转机械设备中[1-2]。目前，人们对滑动轴承转子系统的油膜和运行稳定性等问题一直作为重点研究，在多数的研究中，人们都是通过对简化的雷诺方程求解获得油膜压力，而简化的雷诺方程忽略了N-S方程中的惯性项和油膜曲率等影响因素[1-2]。在文献[3]中作者对抽油机中滑动轴承油膜压力进行了研究，分析了不同偏心率和入口压强对滑动轴承油膜压力的影响，但是并未对所建立的模型进行验证。文献[4]中对不同结构形式的滑动轴承建立模型，通过有限体积法求解不同偏心率和不同转速下的油膜压力分布，并对仿真结果与先前文献中的计算结果进行比较，具有一定的误差。文献[5]建立了滑动轴承有限元模型，通过FLUENT对不同转速下油膜力场进行仿真，而文中只是将仿真结果与弗洛贝格理论与实验曲线趋势进行比较，由于所建模型具有一定的差异，所以结果并不能完全验证结果的准确性。文献[6]建立了椭圆形轴瓦数值分析模型，分析了椭圆度、偏心率等参数对滑动轴承动特性的影响。但是文中只是计算了几种参数对滑动轴承特性影响，而并未对结果的准确性进行验证。文献[7]通过建立椭圆轴承动力特性求解模型，分析椭圆度、偏心率以及转速等因素对滑动轴承油膜压力和承载力的影响，并且比较了椭圆轴承与圆轴承动力特性的差异，文中通过仿真结果与其他研究人员研究结果对比对求解模型进行了验证，其准确性难以保证。

以实验室径向滑动轴承为研究对象，建立相同参数的有限元模型和数值模型，通过有限体积法和有限差分法求解油膜压力场的分布，分析不同偏心率、转速下滑动轴承压力的分布以及油膜压力极值的变化，并对两种方法的计算结果进行对比分析，另外还通过数值模型计算了稳定性临界曲线，并与先前研究人员的结果以及实验结果进行对比，验证模型的准确性和可行性，为进一步研究其他特性提供依。

2 基于有限体积法的模拟计算

2.1 物理模型的建立

以实验室径向滑动轴承为研究对象，建立如图1所示的简化结构模型，结构尺寸参数，如表1所示。

图1 滑动轴承结构模型Fig.1 Hydrodynamic Bearings Structure Model

图中：L—轴承长度；Db—轴承直径；Dj—轴颈直径；c—轴承半径间隙；d—进油孔直径，进油口在轴瓦中心面正上方，O—轴瓦中心线；O1—偏移中心线；e—偏心距；h—油膜厚度。

表1 滑动轴承结构参数Tab.1 Geometrical Parameters of Hydrodynamic Bearings

2.2 油膜网格划分

通过前处理软件对滑动轴承油膜建立三维模型并进行网格划分，如图2所示。

图2 滑动轴承网格结构图Fig.2 Finite Element Grid of Hydrodynamic Bearing

由于轴颈存在一定的偏心，所以滑动轴承模型的油膜厚度较薄，该处网格的精度对最终的计算结果影响较大。若网格划分较大，会使网格产生较大的扭曲率及产生最小体积为负的情况，使网格精度较差，进而影响计算；若网格划分较小，会增加网格的数量，增大分析误差，同时增加计算时间。为提高网格划分精度、计算效率和后期计算数据的提取，将油膜划分为四个区域，分别对各个区域网格划分，根据网格所处位置的不同采用合适的网格密度，并选择轴承间隙径向5层网格，使网格精度更高，最终划分完成的网格，还应能保证能量的传递。考虑网格的精度和计算效率，所划分的网格总数量控制在（50～60）万之间。

2.3 模型假设与边界条件

滑动轴承内的流体为不可压缩流体；润滑油的惯性力忽略不计；润滑油各向同性，黏度为常数；润滑油由三个进油口进入轴承间隙，经计算Re小于2300，油膜间的流动为层流；油膜与轴颈、轴瓦之间没有相对滑移。进油口和出油口均设置为压力边界，压力大小与大气压相等，设置为101325Pa，轴颈表面设置为旋转面，其他设置为固定壁面边界。

2.4 模拟参数设置

通过Gambit网格划分，将网格模型导入到Fluent后，设定润滑油黏度u=0.017Pa·s，润滑油密度为891kg/m3，计算时忽略黏度随温度的变化。选择压力基隐式求解器，离散格式为二阶迎风，算法为SIMPLEC，计算过程中的方程残差为1×10-4，并在计算式对进出口流体质量监控，确定计算收敛。

3 基于有限差分法的数值计算

3.1 滑动轴承油膜厚度

普通滑动轴承油膜厚度可表示为[8]：

由于轴承径向间隙c相对于轴承半径rj为微小量，故式（1）可简化为：

式中：c—轴承径向间隙，c=（Db-Dj）/2；ε—量纲一偏心率，ε=e/c；

Db—轴承直径；Dj—轴颈直径；θ—位置角。

3.2 油膜力求解

非线性油膜压力主要是通过有限差分法对Reynolds直接求解获得。当滑动轴承运行环境满足一定的条件时，Reynolds方程可简化为：

式中：x、z—圆周方向和轴向坐标；p—油膜压力；ρ—油膜密度；

μ—流体动力学黏度；Ut—轴颈t时刻表面切向速度。

在转子稳定运行的过程中，润滑油可认为不可压缩，故润滑油的动力学黏度和密度为固定值，因此，雷诺方程可做进一步简化[9]：

通过采用有限差分法对油膜进行划分网格，并计算任意一个节点处的油膜厚度，将所得油膜厚度带入上式，可计算出每个节点处的油膜压力P。通过数值积分计算出油膜力在U，V方向的分力 Fbu、Fbv：

根据坐标转换可得水平方向与竖直方向油膜力公式：

3.3 系统稳定性临界曲线

通过求解的油膜力，可以求出轴颈任意时刻的运动轨迹，并通过观察不同转速下的轴心轨迹从而确定轴承的临界转速。确定系统的临界转速后，可以绘制一种无量纲的临界曲线图，此稳定性图的横坐标为无量纲的偏心率，纵坐标为无量纲的运行参数，无量纲性的运行参数可表示为[8]：

式中：F—作用于转子上的载荷；棕—转子角速度。

3.4 试验方案

为验证理论的准确性，设计了一种对称刚性转子系统试验台，如图3所示。试验台包括驱动系统、油路系统、信号测试模块和试验平台。滑动轴承转子是由变速电机驱动，通过联轴器联接。在转子x、y方向各有一个非接触式位移传感器，可以测量转子的轴心轨迹。通过位移传感器将采集的位移信号传输到lab VIEW数据记录系统，经过数据分析和处理，可以输出各种类型的试验结果（转子转速、轴心轨迹、振动频谱图等），通过对试验结果分析可获得转子的临界转速，如图4所示。当改变系统的某一参数（转子质量、流体动力学黏度）时，可以测得不同的临界转速，通过计算构成一条无量纲的稳定性曲线[11]。

图3 实验装置Fig.3 Experimental Rig

图4 数据记录系统屏幕Fig.4 View of Data-Logging System

4 结果分析与方法验证

通过运用有限体积法和有限差分法求解了普通滑动轴承油膜压力场，并对其运算结果进行比较分析。

4.1 不同参数下油膜压力分布比较

在轴承结构参数相同，转子的旋转速度为3824r/min时，偏心率着=0.4，轴承油膜压力模拟分布，如图5所示。

图5 油膜压力云图Fig.5 Nephogram of Oil Film Pressure

由图5可知，油膜压力的分布存在2个比较明显的压力区，其中一个为正压区，一个为负压区，且这两个区域分布于最小油膜厚度的两侧，这是由于滑动轴承旋转过程中，在最小油膜厚度两侧，会出现收敛区和发散区。在收敛区，随着间隙不断减小，油膜向内流动形成锲形，使得产生正油膜压力，并在到达最小油膜厚度之前形成最大正压力；在发散区，随着间隙的不断增大，油膜向外流动，压力减小，出现负压区，而在负压区，油膜破裂[10]。

不同参数状态下，轴承Z（轴线）方向中心界面的压力分布，并对模拟结果和理论结果进行了比较，如图6所示。

图6 模拟结果与理论结果油膜压力分布比较Fig.6 Comparison of Simulation Results and Theoretical Results for Oil Film Pressure

由图6可知，给出了模拟结果与理论结果的油膜压力的分布曲线，在不同偏心率和不同转速下，两种计算结果压力分布曲线趋势基本相同，而大小有所不同；在出现油膜压力最大值区域，随转速的增加，油膜正压最大值逐渐增加。在油膜压力开始形成至出现最大油膜压力的过程中，两种计算结果的变化趋势基本一致，而在油膜压力减小的过程中，两种方法的计算结果变化趋势有一定的差异。通过对模拟计算结果分析，在进油孔附近油膜压力分布出现波动，这是由于油膜压力受到进油孔的供油压力的影响。

4.2 不同参数下油膜压力极值比较

模拟结果和理论结果在不同偏心率下油膜压力极值的变化曲线，如图7（a）所示。模拟结果和理论结果在不同转速下油膜压力极值的变化曲线，如图7（b）所示。当偏心率一定时，随着轴颈旋转速度的逐渐增加，油膜压力极值也逐渐增大，增大趋势基本呈线性增加；当转速一定时，随着偏心率的增大，油膜压力极值逐渐增加，且增加速率是逐渐增大的，这是由于偏心率的增大会使轴承的挤压效应增加[3]。由于边界条件、进油压力等差异，模拟结果均稍微高于理论计算的结果。

图7 油膜压力极值的比较Fig.7 Comparison of Oil Film-Pressure Extremes

4.3 无量纲的稳定性曲线比较

通过建立的理论模型计算的一种无量纲的稳定性临界曲线图，如图8所示。纵坐标为无量纲的运行参数op表示稳定性临界转速，其值由式（9）确定，横坐标表示无量纲的偏心率。位于曲线上方的区域为稳定性区域，位于曲线下方的区域为不稳定性区域，而位于曲线上的点为系统标准稳定性临界转速[11]。作者把这里的方法计算的稳定性曲线图与先前研究人员的曲线图以及实验结果进行对比，可见，曲线变化趋势基本一致，证明建立的数值模拟和理论方法在研究滑动轴承特性方面具有可行性。

图8 无量纲稳定性曲线对比Fig.8 Comparison of Stability Map

5 结论

以实验室径向滑动轴承为研究对象，建立了滑动轴承数值计算模型和有限元模型，计算了偏心率和转速对油膜特性的影响以及稳定性临界曲线，得如下结论：

（1）在相同工况下，通过两种方法求解出油膜压力，并进行对比分析，得出油膜压力分布曲线变化趋势基本相同，而大小不同。随着转速的升高，油膜压力极值逐渐增大，模拟仿真结果和理论结果基本呈线性增长；随着偏心率的增大，油膜压力极值增大，模拟仿真结果和理论结果同步增加，且增长速率也逐渐增加。

（2）通过建立的数值计算模型，求解出油膜压力，并计算了滑动轴承转子系统的稳定性临界曲线，将计算曲线与先前研究人员的理论曲线以及实验结果对比，曲线变化趋势基本一致，表明所建模型具有一定的可行性和合理性，为下一步研究形位误差对滑动轴承的影响提供依据。