杨欣，顾根香，孙思聪，周建明，李赛力

1 上海船用柴油机研究所，上海 200090

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

0 引 言

滑动轴承被广泛应用于舰船等各类机械系统中，作为主要的运动副和激励源，它是功率损耗和振动的主要来源之一，动力学及摩擦学特性对其性能及可靠性影响很大[1-2]。因此，对滑动轴承摩擦学和动力学特性进行试验与计算分析具有重要意义。滑动轴承摩擦学和动力学性能研究的3个最重要因素是运动部件的动力学特性、接触面间的摩擦学特性以及其试验验证，但基于有限元法的建模往往会产生较多的自由度，而接触面的摩擦特性也会引入非线性从而导致计算资源消耗的进一步增加。此外，针对发动机内部的轴承的试验条件比较苛刻[3-4]。

常规的有限元离散化法由于会产生过多的自由度而导致求解十分困难，为了在保持计算精度的同时考虑油膜的非线性，近年来发展了动态子结构法，即在建立复杂结构系统力学模型时，从量级上大幅缩减整体结构的自由度，从而保证计算精度[5-7]。针对轴颈和轴瓦间的油膜作用，基于Hertz接触模型提供了一种简化的非线性计算方法[8]。Greenwood等[9]和Patir等[10]分别对弹性粗糙表面的粗糙接触机理和粗糙表面间的油膜动压润滑机理进行研究，得到了Greenwood/Tripp粗糙接触模型和基于平均流量系数的弹性流体动力润滑（EHD）模型，奠定了弹性流体动力润滑理论基础。Bukovnik等[11]等比较了不同类型滑动轴承仿真模型和油液润滑条件的影响。卢伯聪等[12]对主轴承热弹性流体动力润滑（TEHD）进行研究，并在TEHD求解中考虑了轴瓦热效应的影响，从理论上更接近轴承实际的工作状态，但相关研究主要集中在理论计算分析层面，缺乏试验数据的支撑。

本文拟采用仿真与试验相结合的方法，首先基于子结构法建立各零部件的体单元，并通过EHD模型考虑连接副的摩擦学特性，进而建立考虑滑动轴承油膜作用的多柔性体主轴承磨损试验台动力学模型；接着，进行动力学和摩擦学的仿真计算并与试验结果进行对比分析，验证计算方法的正确性，然后在此基础上揭示轴承转子系统的动力学与摩擦学耦合机理；最后，对比基于非线性弹簧模型和EHD模型的计算精度，进而得到主轴承热弹性流体润滑计算的精确建模方法。

1 基本理论

1.1 动态子结构缩减理论

采用有限元方法求解往往会产生规模巨大的离散化模型，加大了求解的难度。子结构法是通过将结构由整体划分为简单的子结构，然后基于特定的算法（静态缩减、古典动态子结构法、约束模态综合法或自由模态综合法等）将原本的大自由度系统缩减为保留主自由度及包含质量、刚度、恢复矩阵的子结构模型。静态缩减往往会使固有频率偏高，且其假设低频段惯性力对副自由度的影响比静力效应小，因此在多体系统动力学中，采用动态缩减进行求解的精度更高，其运动方程为[12]：

式中：M，K分别为模型的质量和刚度矩阵；x，f分别为位移和力矢量；下标m，s分别为主自由度和副自由度的坐标。

由式（1）中的第2个方程式，可以得到

其解为：

取其第i阶特征值λj2的值Pj进行计算，式（3）变为

式中：I为单位矩阵；Ts为坐标转换矩阵。将式（4）代入式(1)，变换后缩减模型的运动方程为

1.2 Greenwood/Tripp粗糙接触模型

随着供油条件变差，当油膜厚度减小至一定程度时，两粗糙表面间部分微凸体直接接触，从而产生较大的摩擦力。针对这种现象，Greenwood等[9]提出了一种考虑表面粗糙形貌名义平面的弹性接触理论，得到了接触变形与表面形貌的相互关系。对于微凸体高度呈标准正态分布的表面，根据Greenwood/Tripp粗糙接触理论，微凸体接触压力pa(h)为

式中：K为弹性因子；E'为等效弹性模量，Pa；Hs为无量纲凸峰间隙；h为润滑油名义油膜厚度，m；为平均凸峰高度，m；σ为两接触表面综合粗糙度，m；ε1，ε2分别为两接触面材料的泊松比；E1，E2为两接触面材料的弹性模量，Pa。

1.3 粗糙接触表面平均雷诺方程

联立流体连续性方程和运动方程（Navier-Stokes方程），假设滑油密度、粘度在油膜厚度方向为常值，可得到简化后的一般形式为雷诺方程。为了考虑流体膜厚度达到与固体表面粗糙度接近的数量级时表面微观形貌对润滑性能的影响，以一般形式的雷诺方程为基础，Patir等[10]引入了压力流量系数、剪切流量系数概念，得到平均雷诺方程为

式中：x，z分别为轴向和径向坐标轴；ϕx，ϕz分别为x，z向压力流量系数；ρ为流体密度，kg/m3；η为机油粘度，mPa·s；h为油膜平均厚度，m；p为油膜压力，Pa；θ为机油填充率；U为轴颈表面线速度，m/s；ϕs为剪切流量系数；t为时间，s。

1.4 EHD数学模型的计算流程

基于动态子结构弹性流体动力润滑（EHD）的计算需联立能量方程、平均雷诺方程、弹性位移方程及载荷平衡方程，并进行求解[8]。其中，能量方程、弹性位移方程采用有限元法求解；平均雷诺方程采用有限差分法求解，且差分网格由轴瓦内表面的缩减模型主自由度节点决定。在计算过程中，采用平均雷诺方程求解得到油膜压力（包含粗糙接触产生的接触压力）后代入弹性位移方程，从而得到各个节点的弹性位移和油膜厚度并迭代至收敛。其计算流程如图1 所示。

图1 EHD仿真计算流程图Fig.1 Flow chart of EHD simulation

2 轴承磨损试验台EHD建模

基于动态子结构缩减理论，利用AVL_Excite Power Unit模块，以轴承磨损试验台为研究对象进行轴承热弹性流体润滑特性分析。轴承磨损试验台主要包含驱动电机、传动轴系、扭矩传感器、被测单元以及液压和供油系统。其中被测单元包含2个支撑轴承座（每个轴承座含2个支撑轴承）、1个被测轴及连杆；驱动电机通过传动轴为被测单元的被测轴提供转速边界条件；液压系统通过液压缸将载荷施加于连杆；供油系统为系统各轴承提供稳定的供油边界条件。轴承磨损试验台的三维结构如图2所示，其基本参数如表1所示。

表1 轴承磨损试验台主要参数Table1 Main parameters of bearing wear test bed

图2 轴承磨损试验台的三维结构Fig.2 Three-dimensional structure of bearing wear test bed

2.1 多体动力学建模

在实际工作中，除被测单元外，其他部件无宏观运动，建模时，只需将被测部分的固定件和运动件视为柔性体，液压载荷以外载荷形式输入，右侧电机的作用通过右侧节点的恒定转速进行考虑即可。试验台包含两侧的支撑座共4个滑动轴承以及被测轴承，因此，所建的动力学模型包含2个轴承座和被测偏心轴以及连杆体单元，其中在偏心轴与各轴承座之间建立2个EHD连接副，在偏心轴与连杆之间建立1个EHD连接副，得到的多体系统动力学模型如图3所示。

图3 轴承磨损试验台动力学模型Fig.3 Dynamic model of bearing wear test bed

2.2 缩减模型的建立

偏心轴和连杆是轴承磨损试验台多体系统动力学模型的主要运动部件，其动力学特性的准确性将直接影响计算结果的精度。为了提高计算效率，采用自由界面子结构法（Craig-Chang法）建立了缩减模型，如表2所示。为了保证缩减模型的动力学特性与有限元模型一致，分别对有限元模型和缩减模型进行了模态计算，其各阶模态频率的对比如表3和表4所示。经对比分析可知，偏心轴、连杆缩减体的各阶模态频率与有限元模型的误差不超过2%，说明缩减模型与有限元模型的动力学特性基本一致。

表2 主要零部件有限元模型与缩减模型Table2 Finite element model and condensed model of main parts

表3 偏心轴有限元模型与缩减模型模态对比Table3 Comparison of modal frequency between finite element models and condensed model of eccentric shaft

表4 连杆有限元模型与缩减模型模态对比Table4 Comparison of modal frequency between finite element model and condensed model of connecting rod

3 试验验证及结果分析

3.1 测试系统

为了验证计算模型的正确性，在轴承磨损试验台布置了相应的测点进行试验数据采集。首先，通过扭矩仪间接测量轴承的摩擦功耗，如图2所示。图4（a）所示为主要测点布置图。在连杆两侧伸出的结构上安装夹角为90°的位移传感器，用于测量轴心轨迹，其原理如图4（b）所示，轴承磨损试验台的测量参数如表5所示。

表5 轴承磨损试验台主要测试参数Table5 Main test parameters of bearing wear test bed

图4 轴承磨损试验台动力学模型Fig.4 Dynamic model of bearing wear test bed

3.2 计算结果验证与分析

3.2.1 结果验证

模型的验证从2个方面进行：1）对比360°曲柄转角内轴承磨损试验台摩擦力扭矩的计算值与试验值，验证其摩擦学特性；2）对比计算与试验所得的轴心轨迹，验证轴承磨损试验台模型转子的动力学特性，如图5所示。

图5 轴承磨损试验台计算结果验证Fig.5 Verification of calculation results of bearing wear test bed

经对比发现，试验测得的扭矩值大于计算值，这是因为试验台在安装过程中存在一定的轴系不对中情况，产生了部分阻力，且自由端储存的滑油也会对轴产生拖拽力矩，不过在摩擦力力矩方面计算值与试验值的误差在10%以内；计算所得与试验测得的轴心轨迹轮廓及位置基本一致，验证了该计算模型在摩擦学及动力学上的正确性。

3.2.2 结果分析

运用本文第2节的建模方法，建立包含5个轴承的轴承磨损试验台EHD多体系统动力学模型，按X方向，从左到右,各轴承号分别为1#～5#。在验证模型精度的基础上，对各轴承的润滑特性进行计算分析。表6对比列出了5个轴承在3 000 r/min工况下的最大油膜压力pmax、最小油膜厚度hmin、最小膜厚比λ、一周期内平均粗糙接触比例γ及摩擦功耗Wm。

表6 各轴承EHD润滑计算结果Table6 EHD calculation results of diffrent bearings of the test bed

由表6可知，1#～5#轴承的润滑状态参数关于3#轴承（被测轴承）对称，被测轴承承受了最大的载荷，且最大油膜压力较其他4个支撑轴承高出一个量级，最小油膜厚度达到了综合粗糙度量级，使得膜厚比偏小，出现了局部粗糙接触。根据Stribeck曲线[13]，两侧支撑轴承最小膜厚比大于10，处于完全液动润滑状态，而被测轴承则处于混合摩擦状态，轴颈与轴瓦间的粗糙接触比例和摩擦功耗呈现相同的规律。被测轴承最小油膜厚度出现的时刻与活塞上的止点相对应，此时，倒置的连杆处于最高位置，下轴瓦与偏心轴表面接触，液压系统对连杆施加的载荷处于峰值。由于支撑轴承的数量为被测轴承的4倍，单个轴承的轴瓦承受的载荷远小于被测轴承下瓦，此时，被测轴承下瓦的油膜压力远大于支撑轴承的油膜压力，且轴颈与轴瓦间出现了部分粗糙接触，产生了粗糙接触压力，进而使轴瓦进入混合润滑状态，粗糙接触比例增高，摩擦功耗增高。

分别提取5个轴承在被测轴承产生最小油膜厚度时刻的油膜总压和粗糙接触压力，其分布图分别如图6和图7所示。由图中不难发现，最外侧固定支撑轴承的油膜压力的分布梯度要小于中间支撑轴瓦的油膜压力分布梯度。另由图7所示粗糙接触压力分布彩图，发现支撑轴瓦无粗糙接触压力，说明即便是处于最恶劣的工况，支撑轴瓦仍可保持热动润滑状态，而被测轴瓦下瓦的两侧则出现了局部粗糙接触。

图6 各轴承油膜总压分布彩图Fig.6 Oil film total pressure distribution of different bearing

图7 各轴承油膜粗糙接触压力分布彩图Fig.7 Oil film asperity contact pressure distribution of different bearings

由上述分析可知，被测的3#轴承其润滑条件最差。为了解其润滑特性，分析不同间隙被测轴承的动力学及摩擦学特性，对比了不同间隙（1～1 000 µm）下，一个周期内的最大油膜压力、最小油膜厚度和平均摩擦功耗。

由图8（a）和图8（b）可以看出，随着半径间隙的增大，轴承最大油膜压力随之升高，但最小油膜厚度在一定范围内是先增大后减小，当半径间隙达20 µm时，最小油膜厚度达到峰值，摩擦功耗最低。当半径间隙非常小时，轴颈与轴瓦间出现微凸体接触，轴承处于边界润滑状态，不易形成足够厚的油膜，而随着间隙的增大，最小油膜厚度逐渐增大，油膜支撑力增强，但当间隙过大时，在端泄的作用下，油膜支撑能力将下降，最小油膜厚度也随之下降。同理，摩擦功耗在半径间隙非常小的时候因边界润滑摩擦系数较大，随着间隙的增大逐渐处于液动润滑状态，摩擦力下降；而当间隙过大时，液动润滑剪切应力会升高，同时由于油膜支撑力下降，局部出现的粗糙接触也会导致摩擦功耗上升。

图8 润滑特性随间隙变化的趋势Fig.8 Variation of lubrication characteristics with clearance

间隙为10和1 000 µm这2种工况下的轴承总摩擦功耗及轴承载荷如图9所示。由图可知，两种工况下的轴承载荷比较接近，但在载荷极值处，即液压系统载荷峰值处，后者较前者力的幅度增大了几百牛，总摩擦耗功也同样从10 µm时的0 W上升到了几百瓦，且摩擦功耗峰值点与载荷极值点对应，说明间隙过大产生的冲击导致局部润滑状态突变，进而使摩擦功耗增加。

图9 一周期内摩擦功耗和载荷Fig.9 Friction power consumption and load in one cycle

3.3 不同建模方式计算结果对比

由上节分析可知，轴承动力学与摩擦学之间存在较强的耦合关系。除上节叙述的通过求解平均雷诺方程进行轴承动力学分析之外，轴承油膜的建模还有其他常用方法，比如简化的非线性弹簧单元，即将轴颈和轴瓦间的高压油膜简化为轴颈与轴瓦间的若干非线性弹簧阻尼单元，求解时，通过激活处于压缩状态的弹簧来迭代计算得到相应的弹簧力，进而得到轴承载荷。为了对比基于EHD模型和非线性弹簧模型建模方法的差别，在相同动力学模型的基础上修改被测轴承连接副的单元类型为非线性弹簧单元并进行计算。将非线性弹簧单元与不同间隙EHD连接副算得的被测试验台轴承动力学、摩擦学结果进行对比，其时域图如图9（b）所示。由图中不难发现，当轴承间隙较小时，采用非线性弹簧模型与EHD模型计算得到的轴承载荷基本一致，但随着轴承间隙的增大，油膜因端泄的作用其支撑力下降；轴颈与轴瓦冲击的增大在载荷时域图中表现为明显的波动。非线性弹簧不存在实际工况中可能出现的冲击现象，采用计算得到的结果在间隙较大的情况下依然是较顺滑的载荷曲线，这与实际工况差异较大。对采用2种建模方法计算得到的轴承载荷进行了傅里叶变换，结果如图10所示。由图可知，这2类模型在低频段重合较好，但采用弹簧单元计算得到的高频分量幅值较EHD单元要低，故对于轴承动力学的求解，使用EHD模型可以提高振动响应的求解精度。

图10 2种建模方式下载荷频域力的对比Fig.10 Comparison of load force between two modeling methods in frequency domain

4 结 论

本文以轴承磨损试验台为研究对象，基于子结构法建立了多体系统动力学模型，然后采用仿真计算与试验相结合的方法研究了滑动轴承摩擦学与动力学的耦合关系，主要得到如下结论：

1） 采用子结构法建立的轴承磨损试验台多体系统动力学模型，通过对动力学和摩擦学分别进行的试验验证，证实了建模方法的正确性。

2） 在试验工况下，轴承磨损试验台的4个支撑轴承处于完全液动润滑状态，被测轴承处于混合润滑状态，这种恶劣的润滑条件主要是因上止点的载荷较高且受力集中所致。

3） 随着间隙的增大，被测轴承的润滑状态先改善并在20 µm间隙时达到最佳，随后恶化；随着间隙的增大，最大压力随之增大，最小油膜厚度先增加后减小，摩擦功耗则是先减小后增大，因此在进行轴瓦设计时，需确定合适的间隙以同时满足润滑性能需求。

4） 当间隙较小时，采用非线性弹簧可以得到较高精度的轴承载荷，但当轴承间隙较大时，在外载的作用下轴颈与轴瓦之间存在冲击，采用EHD模型的求解精度更高。