基于多物理场耦合的可倾瓦滑动轴承性能优化分析*

胡 杨 孟永钢

(1.上海大学机电工程与自动化学院 上海 200072；2.清华大学摩擦学国家重点实验室 北京 100084)

可倾瓦径向滑动轴承由于其良好的稳定性，在大型高速涡轮机械中得到了广泛的应用。当轴径表面线速度较低时，可倾瓦滑动轴承运行在层流状态，可以采用DOWSON[1]提出的广义雷诺方程来预测轴承性能。

然而，在大型涡轮机械中，由于转子直径大、转速高，轴颈表面线速度可达数百米每秒，此时，可倾瓦滑动轴承运行在湍流状态。NG和PAN[2]提出了剪切驱动理论来模拟薄膜流动中的湍流现象。基于NG和PAN[2]的线性化湍流理论，TANIGUCHI等[3]对大型可倾瓦滑动轴承进行了热流体动力学分析，瓦面温度、摩擦损失和偏心率的实验结果与仿真结果吻合较好。

当负载较大时，施加在瓦块表面上的压力使瓦块产生弹性变形；另一方面，温升也会使瓦块产生较大的热变形。轴瓦变形会引起油膜分布的变化，从而影响润滑性能。因此，在重载条件下，有必要考虑轴瓦变形对润滑性能的影响。ETTLES[4]在预测小型轴承的润滑性能时，采用了一维梁弯曲理论来评估轴瓦变形。虽然一维[4]和二维瓦块变形模型[5]具有较高的计算效率，但是当可倾瓦滑动轴承的尺寸较大时，为了确保精度，采用三维瓦块变形模型[6]是必要的。SANO等[7]通过有限元方法评估瓦块变形，将TANIGUCHI等[3]的热流体动力润滑模型扩展为热弹流体动力润滑模型，并基于所构建模型，对直径890 mm的直接润滑的两瓦可倾瓦滑动轴承进行了研究。ARIHARA等[8]构建了热弹流体动力润滑模型，对高速重载条件下的可倾瓦滑动轴承的静态性能进行了研究，并与实验结果进行了对比分析。MERMERAS等[9]描述了一种先进的建模方法，并对新型900 mm三瓦可倾瓦径向滑动轴承进行优化设计和验证。HAGEMANN等[10]在不同供油流量条件下对定向润滑五瓦可倾瓦滑动轴承进行了实验研究，通过与理论结果的对比发现速度对从富油到贫油的转变影响很大。SAN ANDRÉS等[11]评估了在润滑油量过高和过低时可倾瓦滑动轴承的稳态和动态性能。

上述可倾瓦滑动轴承的热弹流体动力润滑模型大多是采用有限差分法或有限元法进行自主编程求解。由于编程代码非开源和不易复制，限制了仿真工具在工业中的广泛应用。使用商业软件可以使研究人员更多地关注物理现象，而不是复杂的数值求解方法。LOHNER等[12]通过结合COMSOL Multi-physics和MATLAB实现了HABCHI[13]的方法，提供了一种求解涂层圆柱线接触的热弹流体动力模型的新思路。受这些工作的启发，本文作者采用了一种结合自主编程和商用软件优点的仿真方法。

本文作者采用COMSOL和MATLAB相结合的方法，建立了考虑湍流的可倾瓦滑动轴承的高效热弹流体动力润滑模型。由于计算流体力学计算量大，且在恶劣工况下易发生失稳，文中所建立的模型仍然是基于雷诺方程。模型中，采用COMSOL中的偏微分方程物理场求解考虑湍流的修正雷诺方程和能量守恒方程以及热传导方程，采用COMSOL中的热应力物理场求解压力和温度引起的瓦块变形，利用MATLAB自编程序来实现子模型之间的耦合。在此基础上，基于所构建模型，研究了轴承瓦块弧长、名义间隙、支承点偏移量和轴向宽度对可倾瓦滑动轴承性能包括油膜厚度、压力、瓦面温度、摩擦功耗等的影响。

1 考虑湍流的可倾瓦滑动轴承热弹流体动力润滑模型

1.1 油膜厚度描述

文中以富油润滑的四瓦可倾瓦滑动轴承作为研究对象，可倾瓦滑动轴承的示意图如图1(a)所示，坐标系如图1(b)所示。油膜厚度表达式如下：

图1 可倾瓦滑动轴承示意与坐标系

h(θ,z)=c+eXcosθ+eYsinθ-m·c·cos(θ-

θp)+δp·(R+tp)·sin(θ-θp)+δd(θ,z)

(1)

式中：c为径向间隙；m为预紧；θp为支承点位置；δp为瓦块倾斜角；R为转子半径；tp为瓦块表面与支承点位置的距离；δd为瓦块变形引起的油膜厚度变化；eX和eY分别为沿X和Y方向轴心位置分量。

1.2 考虑湍流的修正雷诺方程

为了考虑湍流效应的影响，采用NG和PAN[2]推导的线性湍流润滑理论。修正雷诺方程公式表达式如下：

(2)

其中，G1、G2和F1的表达式如下：

(3a)

(3b)

(3c)

(4)

式中：c为径向间隙；N为转速(r/min)；ψ为间隙比。量纲一化坐标系(θ,η,ζ)定义如下：

θ=x/R,η=y/h,ζ=z/(L/2)

(5)

ξ1、ξ2、ξ3和ξ4定义如下：

(6a)

(6b)

(6c)

(6d)

式中：fc(y)和gc(y)为湍流函数，具体计算可参考文献[3]。

1.3 三维能量方程

湍流状态下，量纲一化的三维能量方程表示如下：

(7)

能量方程中量纲一化因子定义如下：

(8)

量纲一化的速度场表达如下：

(9a)

(9b)

(9c)

1.4 热传导方程

量纲一化柱坐标下的瓦块热传导方程表达式如下：

(10)

1.5 瓦块变形

当轴承负载很高时，需要考虑瓦块的变形。瓦块的变形包括两部分，一部分是由于瓦面上油膜压力造成，另一部分是由于瓦块内温度分布不一致造成。瓦块变形量ε的表达式[14]如下:

ε=D-1σ+αT(Tp-Tref)

(11)

式中：σ为应力张量；D是各向同性材料假设下的弹性刚度矩阵；αT为热膨胀张量；Tp为瓦块温度分布;Tref为固体材料热膨胀的参考温度。

1.6 边界条件

热边界条件简述如下：

(1)在流体和瓦块的接触面上，假设热通量是连续的，并且流体和瓦块接触面的温度是相同的。

(2)在流体与轴的接触面上，假设轴表面温度在圆周方向上是相等的，流体与轴之间的整体热交换为0。

(3)在流体入口边界处，温度设定为入口混合温度。根据MITSUI等[15]的研究，使用SUH和PALAZZOLO[16]提出的改进入口混合温度模型来计算瓦块入口的润滑油混合温度，并假设热油携带系数为0.8。

(4)在流体出口边界和轴向两端边界，假设热通量为0。

(5)在其他瓦块表面，假设与周围环境进行热对流。如文献[3]所述，背面的热对流系数为350 W/(m2·K)，其他表面的热对流系数为115 W/(m2·K)。

2 可倾瓦滑动轴承热弹流体动力润滑模型的计算方法

对于某一瓦块内压力场、温度场和瓦块变形场的计算通过COMSOL和MATLAB联合仿真实现。利用COMSOL分别建立3个独立的有限元模型，分别为压力子模型、温度子模型和变形子模型。压力子模型采用COMSOL中偏微分方程物理场模块进行构建，用于求解考虑湍流的修正雷诺方程，从而获得油膜压力。温度子模型也采用COMSOL中偏微分方程物理场模块进行构建，用于求解能量方程和瓦块热传导方程，从而获得流体和瓦块温度场分布。而变形子模型采用COMSOL中热应力物理场模块进行构建，用于获得瓦块变形分布。COMSOL与MATLAB之间的数据传递和交换通过COMSOL的livelink for matlab函数实现[17-18]。在MATLAB中计算偏微分方程的系数，将其写入文本文件，然后在COMSOL中通过插值函数读取。利用COMSOL计算出的油膜压力、瓦块温度和瓦块变形量，通过mpheval、mphinterp、mphglobal等函数输出到MATLAB程序中。此外，COMSOL子模型和MATLAB程序中的网格划分是不同的。在MATLAB中，网格为16×16×10。而在COMSOL中，3种有限元模型有各自的网格划分，如图2所示。此外，压力、温度和瓦块变形的收敛准则采用相对误差为10-3。而力矩平衡和负载平衡的相对误差为10-2。

图2 COMSOL中构建的子模型的网格划分

图3所示为考虑湍流的可倾瓦滑动轴承热流体动力润滑模型总体计算流程。计算开始前，先假定轴心位置和各个瓦块的倾斜角。各个瓦块的力矩平衡通过牛顿-拉斐逊法实现，从而得到各个瓦块的倾斜角。对于无负载的上瓦块可能不存在力矩平衡点。如果没有平衡位置，则给出一个临时倾斜角度(该倾斜角可使力矩的绝对值降至最小)。在给定负载下，轴心位置同样采用牛顿-拉斐逊法迭代获得，直到流体承载力和负载平衡。表1所示为可倾瓦径向滑动轴承结构参数，表2所示为可倾瓦滑动轴承材料参数，表3所示为润滑油参数。

图3 考虑湍流的可倾瓦滑动轴承热弹流体动力润滑模型总体计算流程

表1 可倾瓦径向滑动轴承结构参数

表2 可倾瓦滑动轴承材料参数

表3 HP-8A润滑油参数

3 结果与分析

为了验证所建立的模型，在忽略瓦块变形的情况下，将文中模型预测结果与TANIGUCHI等[3]的预测结果进行比较。仿真中，转速为3 000 r/min，载荷为180 kN，施加在2个瓦块之间。采用直径479 mm、长度300 mm的四瓦可倾瓦滑动轴承。瓦块弧长80°，径向间隙0.612 mm，瓦块厚度121 mm。支承点偏移量50%，预紧为0。润滑油型号为ISO VG32。图4所示为文中模型和TANIGUCHI等[3]模型预测的油膜压力、瓦块表面温度和油膜厚度结果比较。可以观察到，文中的预测结果与TANIGUCHI等[3]的仿真结果吻合较好。

图4 文中模型结果和 TANIGUCHI等[3]仿真结果的对比

与以往自主编程实现的方法相比，采用COMSOL与MATLAB联合仿真的方法可以直观、方便地显示仿真结果，这有助于理解可倾瓦径向滑动轴承运行中的物理现象。以下瓦块2为例，图5所示为可倾瓦滑动轴承下瓦块2 的压力场、温度场和变形场的三维图。其中，载荷为4 927 N，转速为6 000 r/min。红色箭头所示为流体流动方向。如图1所示，载荷施加在下瓦块2和3之间,载荷由下瓦块2和3承受。下瓦块2表面的油膜压力分布如图5(a)所示。由于摩擦产生热量，出口油温高于进口油温，如图5(b)和5(c)所示。靠近瓦块中部的径向变形增大了油膜厚度，而边缘处的径向变形对油膜厚度的影响相反，如图5(d)所示。

图5 可倾瓦滑动轴承瓦块2 的压力场、温度场和变形场的三维图

对于文中研究的四瓦可倾瓦滑动轴承，由于载荷施加在两块下瓦之间，载荷由两块下瓦承受，下瓦表面油膜压力大，而上瓦表面的油膜压力很小。并且，下瓦的瓦温远高于上瓦的瓦温。因此，下文只分析轴承参数变化对下瓦的影响。仿真中，负载为20 000 N，转速为10 000 r/min。为了研究轴承结构参数对可倾瓦滑动轴承性能的影响，将采用控制变量法进行研究。轴承参数包括瓦块弧长、径向间隙、支承点偏移量和轴向宽度。

3.1 瓦块弧长对可倾瓦滑动轴承性能的影响

图6所示为瓦块弧长对可倾瓦滑动轴承性能的影响。随着瓦块弧长的增加，承载面积变大，面压(F/(L×D))降低，从而油膜最大压力会降低，如图6(a)所示。同时，由于面压降低，偏心率降低，油膜厚度增加；进一步地，由于油膜增厚，油膜局部剪切程度降低，瓦面最大温度降低，如图6(b)所示。图6(c)所示为油膜厚度随瓦块弧长的变化趋势，可见油膜厚度随瓦块弧长的增加而增加。而由于偏心率降低，轴心位置会上浮，如图6(d)所示。

图6 瓦块弧长对可倾瓦滑动轴承性能的影响

表4所示为瓦块弧长对可倾瓦滑动轴承摩擦功耗的影响。

表4 瓦块弧长对可倾瓦滑动轴承摩擦功耗的影响

由表4可以看出，随着瓦块弧长的增加，承载面积增加，虽然由于膜厚增加，油膜局部剪切变弱，但可倾瓦滑动轴承的摩擦功耗也会增加。此外，需要注意的是，由于载荷施加在2个下瓦之间，两块下瓦上的压力、温度、油膜厚度分布基本相同。

3.2 径向间隙对可倾瓦滑动轴承性能的影响

图7所示为径向间隙对可倾瓦滑动轴承性能的影响。随着径向间隙的增大，由于负载不变，油膜压力变化很小，如图7(a)所示。此外，由于径向间隙增大，导致偏心率增加，油膜厚度增加，油膜局部剪切变弱，瓦面最大温度降低，如图7(b)所示。油膜厚度随径向间隙的变化趋势如图7(c)所示，可见油膜厚度随径向间隙的增加而增加。由于偏心率的增加，轴心位置相应下沉，如图7(d)所示。表5所示为径向间隙对可倾瓦滑动轴承摩擦功耗的影响。随着轴向间隙的增加，油膜厚度增加，油膜剪切应力降低，摩擦功耗会略有降低。

图7 径向间隙对可倾瓦滑动轴承性能的影响

表5 径向间隙对可倾瓦滑动轴承摩擦功耗的影响

3.3 支承点偏移量对轴承性能的影响

图8所示为支承点偏移量对可倾瓦滑动轴承性能的影响。从图8(a)可以看出，随着支承点偏移量的增加，瓦块2和瓦块3的油膜压力分布差异性越来越大。随着支承点偏移量的增加，瓦块2的油膜压力增大，而瓦块3的油膜压力降低，如图8(a)所示。与此同时，瓦块2的油膜厚度减小，油膜局部剪切变强，瓦块表面最大温度增大，相反，瓦块3的油膜厚度减小，油膜局部剪切变弱，瓦块表面最大温度降低，进而瓦块2的表面温度大于瓦块3的表面温度，如图8(b)所示。随着支承点偏移量的增大，瓦块2 和瓦块3的油膜厚度的变化趋势相反，如图8(c)所示。而对于轴心位置有向右上方移动的趋势，如图8(d)所示。表6所示为支承点偏移量对可倾瓦滑动轴承摩擦功耗的影响。可以看出，随着支承点偏移量的增加，承载能力会有所提高，可倾瓦滑动轴承的摩擦功耗略微增加。

图8 支承点偏移量对可倾瓦滑动轴承性能的影响

表6 支承点偏移量对可倾瓦滑动轴承摩擦功耗的影响

3.4 轴向宽度对轴承性能的影响

图9所示为轴向宽度对可倾瓦滑动轴承性能的影响。可以看出，随着瓦块轴向宽度的增加，承载面积增加，承载能力增大,而面压降低，油膜最大压力降低，如图9(a)所示。面压降低会导致偏心率减小，从而油膜变厚。由于油膜增厚，油膜局部剪切程度降低，瓦面最大温度会相应降低，如图9(b)所示。油膜厚度随着轴向宽度的变化趋势如图9(c)所示，可见油膜厚度随着轴向宽度的增大而减小。随着偏心率降低，轴心位置有上浮的趋势，如图9(d)所示。表7所示为轴向宽度对可倾瓦滑动轴承摩擦功耗的影响。一方面，随着轴向宽度的增加，接触面积增加；另一方面，由于膜厚增加，油膜局部剪切变弱。相比之下，接触面积增大的影响要比油膜厚度局部剪切变弱的影响显著，因此，可倾瓦滑动轴承的摩擦功耗会大幅增加。

图9 轴向宽度对可倾瓦滑动轴承性能的影响

表7 轴向宽度对可倾瓦滑动轴承摩擦功耗的影响

4 结论

结合商用软件和自主编程的优点，采用COMSOL和MATLAB联合仿真的方法构建了考虑湍流的可倾瓦滑动轴承热弹流体动力润滑模型。通过与文献对比，验证了模型的准确性。基于所构建模型，研究了轴承参数包括瓦块弧长、径向间隙、支承点偏移量和宽度对可倾瓦滑动轴承性能的影响。主要结论如下：

(1)随着瓦块弧长的增加，油膜压力降低，油膜厚度增加，瓦面温度降低，而轴心位置略有上浮，可倾瓦滑动轴承的摩擦功耗也会增加。

(2)随着径向间隙的增大，油膜压力变化很小，瓦面温度降低，油膜厚度增加，轴心位置相应下沉，摩擦功耗会略有降低。

(3)随着支承点偏移量的增加，瓦块2的油膜压力、温度增大，而油膜厚度降低；而瓦块3的油膜压力，温度降低，而油膜厚度增加。轴心位置有向右上方移动的趋势，承载能力、摩擦功耗略微增加。

(4)随着瓦块轴向宽度的增加，油膜压力降低，瓦面温度降低，油膜厚度增加，轴心位置有上浮的趋势，摩擦功耗大幅增加。

对于可倾瓦滑动轴承，瓦面最大温度和摩擦功耗是2个重要的指标。瓦面温度关系到轴承安全运行，而摩擦功耗会影响节能减排。因此，需要对轴承参数进行优化选择。文中的研究为可倾瓦滑动轴承的优化设计提供有价值的参考。