基于格心型有限体积法的滑动轴承液膜压力分析

杨国栋， 曹贻鹏， 明平剑， 张文平

(哈尔滨工程大学 动力与能源工程学院，哈尔滨 150001)

滑动轴承由于形式简单，接触面积和承载能力大，广泛应用于工业生产和生活领域中。当前，滑动轴承动压润滑现象的求解方法主要包括有限差分法[1]、有限元法[2]及CFD软件仿真[3-4]等。有限差分方法(Finite Differential Method，FDM)是计算机数值模拟最早采用的方法，发展也较为成熟，但是FDM对求解域的几何特征要求较高，其结构化网格不易处理复杂的区域，同时在网格数较少时，较难满足积分守恒；有限元法(Finite Element Method， FEM)主要应用于结构力学领域，近些年开始应用于流体力学领域，能够采用非结构化网格处理几何特征复杂的区域，但是FEM需要求解刚度矩阵，当单元数较多时，会占用较大的计算机内存，进而使求解速度下降[5]；一些学者也采用CFD软件仿真的方法来研究轴承的润滑性能，通过求解N-S方程来得到轴承的压力分布，这种方法更为准确，但是求解较复杂，而且油膜厚度是微米级，厚度方向的网格划分较为困难[6]。

李强等[7]基于FVM(Finite Volume Method)研究了考虑气穴现象的滑动轴承油膜压力特性，但是其方程推导是基于结构化网格进行的，这种方法本质上有限差分法一致，没有体现出来有限体积法对复杂网格的适应性。

格心型有限体积法(Cell-Centered Finite Volume Method， CCFVM)是基于CFD技术中的FVM的离散方法建立起来的[8]，所有变量定义在单元中心，要求因变量的积分守恒对任意一组控制体积都得到满足，进而对整个计算区域也满足守恒；该方法对复杂的求解域有很好的适应性，当网格数较多时，处理效率也比较高，因此该方法广泛应用于流场数值计算和力学相关领域，而在润滑领域应用较少。

本文基于CCFVM离散雷诺方程，推导适用于求解非结构化网格的方程离散形式，将计算域的网格要求由结构化网格扩展到非结构化网格，可用于复杂几何表面之间的润滑特性研究，并对文献[9]中的算例进行了验证，在此基础上，以水润滑橡胶轴承为例，分析了长径比、偏心率对其润滑特性的影响。

1 理论方程

1.1 雷诺方程

雷诺方程是由简化的纳维-斯托克斯方程和连续方程推导而得的，常用来求解转子机构的动静压润滑问题，适用于二维、稳态、不可压流体动压润滑计算的雷诺方程形式为[10]

(1)

式中：x为圆周方向坐标；y为轴承宽度方向；h为油膜厚度；η为滑油黏度；U为轴颈转速；P为油膜压力。

1.2 液膜厚度方程

滑动轴承液膜厚度h的计算公式为

h=c+ecosφ

(2)

(3)

式中：c为轴承间隙；e为轴承偏心距；φ为轴承圆周方向角坐标；R为轴承半径。

2 数值方法分析

2.1 雷诺方程的矢量表达形式

雷诺方程式(1)可以用如下形式表示

(4)

对式(4)，在空间上对控制体V进行积分，有

(5)

2.2 扩散项的离散[11]

根据高斯散度定理，将体积分转化为面积分并离散可得

(6)

式中：nf为当前单元周围的单元个数；Af为面积矢量Af的模。

图1 CCFVM控制体示意图Fig.1 Sketch map of the control volume of CCFVM

图1中：ef为两个单元中心的距离矢量；d=ef/ef为单位距离矢量；n为单元面的单位外法线矢量。

由于

(7)

将式(7)代入式(6)中可得

(8)

单元面上的变量梯度由当前单元中心和相邻单元中心的变量梯度线性插值得到。

2.3 源项的处理

式(5)中的源项在控制体内积分可得

(9)

2.4 雷诺方程的离散形式

由式(8)和式(9)可得式(1)的离散形式为

(10)

2.5 边界条件

对于径向滑动轴承，其求解域可展开为矩形，采用雷诺边界条件

(11)

式中：x′为待定的出口边界；B为轴承的宽度。

边界上变量直接赋值为Dirichlet边界条件，边界上变量梯度直接赋值称之为Neumann边界条件。这两类边界条件对式(10)中的贡献为

式中：nD为位于Dirichlet边界SD上所有单元面的个数；nN为位于Neumann边界SN上所有单元面的个数； 下标Bf表示边界面。

2.6 收敛条件

CCFVM算法的收敛条件为

(12)

式中：ξ=1×10-6。

3 算例验证与分析

本文的计算程序是在哈尔滨工程大学动力装置工程技术研究所自主开发的通用输运方程求解器GTEA软件的基础上开发的，采用Fortran语言编程，网格在Gambit 2.4.6中进行划分。

3.1 验证算例

为了验证CCFVM求解二维稳态润滑问题的适用性和准确性，本文将通过Sun等研究中的算例来验证。

表1 径向滑动轴承基本参数

采用非结构化网格划分求解域，共1 556个节点，2 946个三角形单元；孙军采用5 400个节点，5 236个单元。本算例中偏位角为0°，轴颈倾角分别为0°，0.004°，0.007°，0.01°。

计算完毕后，以Sun等论文中的压力结果(图2～图5中左侧的图)为参考，进行以下对比：

由图2～图5和表2可知，当前算法的液膜压力分布同参考值趋势一致，最大压力值基本相当，验证了CCFVM算法在求解雷诺方程及滑动轴承润滑特性方面的适用性和准确性。

图2 倾斜角为0°Fig.2 Slope angle=0°

图3 倾斜角为0.004°Fig.3 Slope angle=0.004°

图4 倾斜角为0.007°Fig.4 Slope angle=0.007°

图5 倾斜角为0.01°Fig.5 Slope angle=0.01°

参数数值倾角/(°)00.0040.0070.01参考值/MPa33.0639.663.58415.35当前值/MPa33.0139.3763.62412.7偏差/%0.150.580.060.72

3.2 实例计算

以水润滑橡胶轴承为例，研究不同工况下的轴承润滑特性。采用CCFVM研究轴承长径比、偏心率的变化对轴承润滑特性的影响。

表3 艉轴承的基本参数

为了研究长径比的变化对轴承压力分布的影响，保持轴承宽度不变，长径比分别为1，2，3，4，计算得到的压力分布如图6所示。

由图6可知，沿轴承圆周方向，液膜压力先增大后减小；随着长径比的增大，高压区沿轴承宽度方向逐渐扩张，压力峰值逐渐下降，说明长径比增大会降低轴承的峰值压力，同时提高轴承的承载力。

为了研究偏心率对液膜压力的影响，保持长径比为1，其他参数不变，偏心率分别取0.1,0.2,0.3，0.4,05,0.6,0.7,0.8,0.9，结果如图7所示。

由图7可知，随着偏心率的增大，液膜压力峰值逐渐增大，偏心率在0.7～0.9变化时，压力增大的更为剧烈；在圆周方向，压力的分布遵循收敛楔的几何变化，先增大后减小。

图6 长径比对液膜压力分布的影响Fig.6 Effects of aspect ratio on film pressure distribution

图7 偏心率对液膜压力的影响Fig.7 Effects of eccentricity ratio on film pressure

4 结 论

(1) 提出了数值求解雷诺方程的CCFVM算法，该算法物理意义更明确，由于采用了非结构化网格处理技术，可以适用结构复杂的润滑区域，同时不会消耗太多存储空间，相比FDM(结构化网格)和FEM(消耗较大存储空间)具有明显的优势。

(2) 通过与文献的结果对比，液膜压力分布和趋势基本吻合，最大液膜压力值接近，验证了该算法的适用性和正确性较好。

(3) 随着倾角的增大，压力峰值逐渐向轴承尾端偏移，同时最大压力值急剧增大；长径比越大，液膜高压区沿轴承宽度方向扩张的越大，同时压力峰值会减小；偏心率越大，轴承压力峰值越大，当偏心率大于0.7时，峰值变化更为明显。