基于弹流耦合算法的船舶艉轴承润滑特性分析*

张 斌 朱汉华

(武汉理工大学能源与动力工程学院 湖北武汉 430063)

大型船舶艉轴承经常处于重载工况，交变载荷产生的摩擦热使润滑油油膜厚度减小，甚至使轴瓦承载区无法形成具有一定承载力的油膜[1]。此时，弹性变形和油膜厚度非常接近，相对误差如此之小的弹性变形势必会对轴承润滑性能产生重要影响。弹流润滑理论的日趋完善，使广大学者渐渐看到了弹性变形的本质和影响。

CONWAY和LEE[2]利用弹性位移方程来研究刚性轴承座中特殊形状轴瓦表面弹性位移分布的变化，并简化了数值计算求解弹性变形的过程。FANTINO和FRENE[3]对表征内燃机连杆轴承的润滑性能参数进行了迭代求解，对比分析了是否考虑弹性变形润滑性能产生的差异，结果表明，计入弹性变形后油膜压力呈现双峰甚至多峰分布，同时最小油膜厚度近似线性减小。何芝仙和桂长林[4]优化了弹性变形的计算方法，虽然结果存在一定误差，却大大提高了计算速度。邹争等人[5]通过数值计算的方法得到艉轴承变形的解析解，分析了安装参数和工况参数对轴承变形的影响，并通过设计正交表来评估不同影响因素对弹性变形的影响程度。马艳艳和崔子伟[6]发现应力偶轴承更适合研究弹性变形，并在雷诺方程中加入应力偶参数，推导出适合分析应力偶轴承变形的雷诺方程。结果表明：相比于普通轴承，应力偶轴承的弹性变形更明显，且最大弹性变形量与应力偶参数成正比。近年来，随着对弹性变形研究的逐步深入，单一变量求解的方法已经不再满足精度要求，研究人员提出了油膜压力、厚度、弹性变形的耦合求解算法，使得弹流润滑理论日趋完善。

本文作者以重载工况下的船舶艉轴承为研究对象，通过有限单元法结合耦合算法求解油膜压力、油膜厚度、弹性变形，进一步分析不同影响因素对艉轴承弹流润滑特性的动态影响趋势，得到一些有利于轴承设计和安装的结论。

1 弹流润滑数学模型

1.1 Reynolds方程

对Navier-Stokes方程(N-S方程)进行简化得到Reynolds方程[7]：

(1)

式中：h为油膜厚度；p为油膜压力；η为润滑油动力黏度；U为轴颈转速；x为周向坐标；y为轴向坐标。

式(1)为研究流体润滑提供了数值计算的依据，其意义在于将油膜压力p、油膜厚度h、润滑油动力黏度η、轴颈转速U等看似没有关联的参数有效联系起来，实现了已知某一参数来求解其他未知参数的目标。文中将基于Reynolds方程进行润滑分析。

1.2 有限单元法求解轴瓦弹性变形

文中建立了考虑轴瓦厚度的三维有限单元模型，将轴瓦整个表面的弹性变形分解到每个节点来求解。对轴瓦模型建立三维坐标系，x、y、z分别表示周向、轴向和径向。所建模型及坐标系如图1所示。

(1)离散和选择一般单元

在对多自由度的计算模型进行离散时通常选用8节点六面体单元，文中轴瓦模型单元用u、v、w分别表示每个节点的3个自由度，下标1～8分别表示六面体的8个节点。如图2所示。

(2)选择位移函数

假设每个节点的坐标为(xi，yi,zi)，i=1,2,…,8

选择位移函数[8]：

u=a1+a2x+a3y+a4z+a5xy+a6xz+a7yz+a8xyz

(2)

各节点的位移分别表示为u1,u2,…,u8，位移矩阵如式(3)所示。

(3)

定义形状函数为N，节点位移量为S，则：

(4)

将3个自由度的位移场u、v、w分别表示为向量形式：

(5)

(3)定义应变位移和应力应变的关系[9]

三维实体模型的应变向量如式(6)所示。

(6)

定义矩阵B(6×24)为

(7)

则应力应变的关系可以表示为

ε=BS

(8)

(4)推导单元刚度矩阵

单元刚度矩阵由矩阵B、D通过转置变换和积分推导求得[10]：

(9)

式中：Ke为单元刚度矩阵；D为弹性矩阵。

(5)组装形成整体刚度矩阵

如果没有引入适当的边界条件对单元刚度矩阵进行自由度限制而直接求解，很难消除多自由度带来的不确定性，进而引发位移结果的多变性。因此，合适的边界条件对于刚度矩阵的求解至关重要[11]。文中引入如下的边界条件：

①三维有限单元模型表面的等效压力载荷来自于轴瓦内表面受到的油膜压力，且油膜压力连续分布在整个轴瓦内表面；

②在组装整体刚度矩阵的过程中，轴瓦外表面受到的约束为刚性约束，模型各节点均无任何自由度的位移[12]。

组装形成的整体刚度矩阵为

(10)

式中：K为整体刚度矩阵；G为单元节点转化矩阵。

(6)求解位移

求解位移，需建立有限元方程进行求解，进而得到整体刚度矩阵中各节点的位移式[13]。文中建立的有限元方程如下：

KSz=Fz

(11)

式中：Sz为整体结构中各节点的位移；Fz为整体结构中各节点的等效载荷。

1.3 弹流润滑耦合算法

弹流润滑耦合算法具体方法为：首先通过有限差分法求解Reynolds方程得到油膜厚度和压力，作为弹流润滑油膜厚度和压力的初值；在初始压力作用下轴瓦会产生弹性变形，计入弹性变形后油膜厚度会改变[14]，此时需要在初始油膜厚度上加入弹性变形，得到新的油膜厚度；油膜厚度的改变又会影响油膜压力，再代入Reynolds方程求解得到新的油膜压力。以上步骤作为第一次循环，通过MATLAB编程进行多次迭代，事先设定收敛精度，每一次循环结束后都校核轴承给定外载荷和油膜承载力的误差[15]，检验误差是否满足收敛精度，这样才能真正实现对润滑性能参数的耦合求解。耦合算法的求解流程如图3所示。

图3 弹流润滑耦合算法流程Fig.3 Flow of elastohydrodynamic lubrication coupling algorithm

2 模型验证

为验证所建弹流润滑数学模型的正确性，选取文献[16]中的全圆轴承参数进行验证，基本参数如表1所示。

表1 轴承及润滑油基本参数Table 1 The basic parameters of bearing and lubricating oil

对比图4和图5可知，虽然文中得到的压力分布图中从峰值压力降为0的过程要比较缓慢，可能是因为在求解轴颈平衡位置的过程中文献[16]是采用割线法，而文中采用非线性方程求根的方法，造成一定误差，但油膜压力沿周向和轴向的变化趋势基本一致，可以证明文中数学模型的正确性。

图4 文中油膜压力分布Fig.4 Oil film pressure distribution in this paper

图5 文献[16]油膜压力分布Fig.5 Oil film pressure distribution in reference[16]

3 动态参数对艉轴承弹流润滑特性的影响

基于建立的船舶艉轴承弹流润滑数学模型，将弹性模量、轴承间隙、长径比3种参数作为影响因素，分析了3种动态参数对艉轴承弹流润滑特性的影响。所选用的艉轴承相关参数如表2所示。

表2 艉轴承及润滑油相关参数Table 2 Parameters of stern bearing and lubricationg oil

3.1 动态弹性模量对艉轴承弹流润滑特性的影响

弹性模量在弹性力学中定义为物体所受的应力和产生应变的比例系数，作为弹性体在弹性变形阶段的重要物理量，可以用来衡量材料产生变形的难易程度，也可以看作材料刚度的体现。图6—13示出了动态弹性模量对艉轴承弹流润滑特性的影响，弹性模量的范围是100～150 GPa。

图6 动态弹性模量影响下弹性变形三维分布Fig.6 Three dimensional distribution of elastic deformationunder the influence of dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

图7 弹性模量对最大弹性变形量的影响Fig.7 Influence of elastic modulus on maximum elastic deformation

图8 动态弹性模量影响下油膜压力三维分布Fig.8 Three dimensional distribution of oil film pressureunder the influence of dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

图9 弹性模量对最大油膜压力的影响Fig.9 Effect of elastic modulus onmaximum oil film pressure

图10 动态弹性模量影响下油膜厚度三维分布Fig.10 Three dimensional distribution of oil film thicknessaffected by dynamic elastic modulus：(a)100 GPa;(b)110 GPa;(c)120 GPa;(d)130 GPa;(e)140 GPa;(f)150 GPa

图11 弹性模量对最小油膜厚度的影响Fig.11 Influence of elastic moduluson minimum oil film thickness

图12 弹性模量对摩擦力的影响Fig.12 Effect of elastic modulus on friction

图13 弹性模量对端泄流量的影响Fig.13 Effect of elastic modulus on end discharge

对比图6和图8可以发现,弹性变形和油膜压力的分布趋势相似，几乎都是在周向中点附近取得最大值，沿周向和轴向均呈现先增后减的趋势。由图10可知，油膜厚度在周向中点附近取得最小值。轴承在周向180°附近油膜压力和油膜厚度会同时出现最大值和最小值，因此该区域润滑环境复杂，需要引起关注。通过观察图7、图9和图11可知，弹性模量增加时，油膜峰值压力增加，最大弹性变形量和最小油膜厚度均减小。原因在于弹性模量增加，轴瓦表面刚性变强，轴瓦不易发生形变，相同外部条件下最大弹性变形量减小；轴承承载能力增加，意味着油膜峰值压力增加；增大油膜压力会导致偏心率增大，最小油膜厚度在周向0～180°的区间内出现，由油膜厚度的公式可知，在此区间内油膜厚度与偏心率成反比，因此最小油膜厚度会减小。

由于黏度的作用，润滑油在轴承间隙流动时会产生摩擦力。如果摩擦力太大，会产生过量的摩擦热和摩擦功耗，影响轴承工作效率。轴承在工作中会从轴承两端泄漏部分油量，称为端泄流量。为保持一种动态平衡，需要不断补充和端泄流量相等的润滑油，所以计算端泄流量的意义在于了解补偿流量。由图12和图13可知，增大弹性模量，材料的摩擦因数增大，因此摩擦力增加；端泄流量与油膜压力和厚度都有关，虽然膜厚减小，但周向两端压力梯度增大的作用更明显，导致端泄流量小幅度增加。

综合弹性模量对艉轴承弹流润滑性能的影响可以看出，在一定区间内增大弹性模量能有效减小轴瓦产生的弹性变形，油膜峰值压力、摩擦力和端泄流量均增加，但最小油膜厚度会近似线性减小。

3.2 动态轴承间隙对艉轴承弹流润滑特性的影响

弹性变形和轴承间隙之间有着不可分割的关系，而间隙改变会直接影响到油膜厚度。此外，间隙太小可能使轴颈和轴瓦直接接触，发生摩擦磨损，导致润滑失效。图14—21示出了艉轴承润滑特性随轴承间隙的动态变化趋势，轴承间隙的范围是0.03～0.08 mm。

图14 动态轴承间隙影响下弹性变形三维分布Fig.14 Three dimensional distribution of elastic deformationunder the influence of dynamic bearing clearance：(a)0.03 mm;(b)0.04 mm;(c)0.05 mm;(d)0.06 mm;(e)0.07 mm;(f)0.08 mm

图15 轴承间隙对最大弹性变形量的影响Fig.15 Influence of bearing clearance onmaximum elastic deformation

图16 动态轴承间隙影响下油膜压力三维分布Fig.16 Three dimensional distribution of oil film pressureunder the influence of dynamic bearing clearance：(a)0.03 mm;(b)0.04 mm;(c)0.05 mm;(d)0.06 mm;(e)0.07 mm;(f)0.08 mm

图17 轴承间隙对最大油膜压力的影响Fig.17 Influence of bearing clearanceon maximum oil film pressure

图19 轴承间隙对最小油膜厚度的影响Fig.19 Influence of bearing clearance on minimum oil film thickness

图20 轴承间隙对摩擦力的影响Fig.20 Influence of bearing clearance on friction

图21 轴承间隙对端泄流量的影响Fig.21 Influence of bearing clearance on end discharge

由图14—19可知，弹性变形、油膜压力沿周向和轴向均近似呈抛物线分布，且基本分布趋势不受轴承间隙变化的影响。轴承间隙改变会直接影响弹性变形量，进而改变油膜压力、厚度，可以看出，当轴承间隙增加时，油膜峰值压力增加，最大弹性变形量和最小油膜厚度则线性减小，与轴承间隙呈现出较强的负相关关系。从图20、21可以看出，摩擦力、端泄流量随轴承间隙增大而增大，但增加幅度不大，均在10%以内。这主要是因为增加轴承间隙和增大材料的泊松比效果相似，会导致轴瓦的竖直变形量增加，因此最大弹性变形量增加。在定偏心距时，由偏心率和轴承间隙的关系可知，间隙增大会引起偏心率减小，而偏心率是相关润滑参数的重要影响因素，其值减小会直接引起油膜峰值压力、摩擦力、端泄流量增加，最小油膜厚度减小，尤其是大偏心率情况下，润滑参数的变化会更明显。

综上所述，轴承间隙增加可以有效减小最大弹性变形，提高油膜峰值压力，但最小油膜厚度会减小，而摩擦力和端泄流量受轴承间隙的影响不大。因此在轴承安装时，需要将轴承间隙控制在合理的区间内，确保轴承处于优良的润滑环境。

3.3 动态长径比对艉轴承弹流润滑特性的影响

长径比作为材料重要的设计参数之一，通常用来描述材料的形状，长径比大的材料为细长状，长径比小的材料为短粗状。图22—29示出了艉轴承弹流润滑特性随长径比的动态变化趋势，长径比的范围是0.5～5。

图22 动态长径比影响下弹性变形三维分布Fig.22 Three dimensional distribution of elastic deformationunder the influence of dynamic aspect ratio：(a)L/d=0.5;(b)L/d=1;(c)L/d=2;(d)L/d=3;(e)L/d=4;(f)L/d=5

图23 长径比对最大弹性变形量的影响Fig.23 Influence of length diameter ratioon maximum elastic deformation

图27 长径比对最小油膜厚度的影响Fig.27 Influence of aspect ratio on minimum oil film thickness

图28 长径比对摩擦力的影响Fig.28 Effect of length diameter ratio on friction

图29 长径比对端泄流量的影响Fig.29 Effect of length diameter ratio on end discharge

从图22、24可看出，弹性变形与油膜压力出现最大值的位置仍在周向中点附近。从图23、25、27可知，当长径比从0.5线性增加至5时，最大弹性变形量近似线性增加，油膜峰值压力则刚好相反，而最小油膜厚度仅增加了0.001 3 mm，说明长径比对最小油膜厚度的影响较小。从图28、29可知，摩擦力和端泄流量都与长径比呈现出一种负相关关系。总体来看，增大长径比，会使轴瓦内表面产生较大的弹性变形，轴承间隙减小，进而降低油膜峰值压力、摩擦力和端泄流量。产生上述变化的原因在于增大长径比，轴瓦受力的比面积减小，相同外载条件下会产生更大的形变，因此最大弹性变形量增加；改变长径比的本质在于改变了轴承形状，大长径比的轴承会使油膜压力分布更均匀，比压不变时，油膜峰值压力减小，最小油膜厚度略有增加。此外，长径比也是偏心率的重要影响因素之一，增大长径比致使偏心率增加，而摩擦力和端泄流量都会随偏心率的增加同步减小。因此，在设计艉轴承长径比时，应综合考虑艉轴承在重载工况下可能产生的弹性变形以及弹性变形对润滑特性的影响。

4 结论

基于弹流润滑耦合算法分析了艉轴承的弹流润滑特性，探讨弹性模量、轴承间隙、长径比3种因素对艉轴承弹流润滑特性的影响，得到以下结论：

(1)弹性变形和油膜压力的三维分布趋势相似，沿周向和轴向都近似呈抛物线分布，最大值出现在周向中点附近。

(2)随弹性模量增加，油膜峰值压力增加，最大弹性变形量和最小厚度均减小，摩擦力和端泄流量同时增加；随轴承间隙增加，油膜峰值压力稳步提高，最大弹性变形量和最小油膜厚度却刚好相反，近似线性减小，而摩擦力、端泄流量受轴承间隙的影响较小；随长径比增加，弹性变形与油膜压力的分布更为平缓，最大弹性变形量近似线性增加，油膜峰值压力减小，最小油膜厚度变化不明显，而摩擦力和端泄流量都与长径比呈现出一种负相关关系。