滚动轴承外圈故障的显式有限元动态仿真分析

时间：2024-07-28

李国超彭炜李勇才高立新张键

1.北京工业大学北京市先进制造重点实验室，北京，100124

2.武汉钢铁（集团）公司，武汉，430081

3.武汉昊海立德科技公司，武汉，430080 4.武汉工程职业技术学院，武汉，430415

0 引言

滚动轴承是现代生产和生活中应用最广泛的零部件之一，其运行状态对设备的正常运转及正常功能的实现起着非常重要的作用。然而轴承也是极易损坏的零件之一，其主要损坏形式包括：疲劳、磨损、腐蚀、电蚀、塑性变形、断裂和开裂等［1］。据统计，由轴承引起的故障占旋转机械设备故障的7%左右［2］，其中90%的故障发生在内圈和外圈上。目前对于轴承故障的研究主要集中在轴承的检测和诊断上，即利用数据采集系统及信号处理技术［3］提取轴承损坏元件的特征频率，以达到对轴承故障进行预测的目的。

用有限元方法对滚动轴承进行仿真分析的有效性已被验证。本文以显式算法［4］为基础，应用ANSYS／LS-DYNA软件对滚动轴承常见的外圈裂纹故障进行了仿真分析。在物理模型的基础上，选取贴合实际的材料模型，并选择合理的摩擦及约束条件对故障模型进行了仿真分析。仿真结果为研究轴承故障机理及寻找轴承故障特征提供了条件和依据。

1 有限元模型的建立

滚动轴承的内部运动相当复杂，不仅存在着滚动体相对内外圈的运动，更有因此而产生的离心力和陀螺力矩，同时还存在着因滚动体与轴承内外圈接触而产生的摩擦力和油膜拖动力［5］。基于滚动轴承上述工作工况，对滚动轴承模型进行如下假设：

（1）因轴承倒角对轴承内部应力分布影响较小，故仿真模型对倒角进行了简化。

（2）未考虑径向游隙和轴向游隙及油膜对轴承运转的影响。

（3）轴承主体均采用线弹性材料。

所分析轴承为实验室轴承试验台6307深沟球轴承，其尺寸参数如表1所示。根据滚动轴承故障特征频率计算方法得到的轴承外圈故障的特征频率为73.1Hz［6］。

表1 滚动轴承6307的实际参数

所建轴承有限元模型主体选用Solid164单元，为增加其旋转自由度，将内圈内表面定义为Shell163单元。故障设置为4mm×1mm的裂纹故障，如图1a所示。网格划分采用扫略、映射和自由划分相结合的方式，划分后的模型如图1b所示。划分的网格共有节点21 030个，单元49 967个。

图1 滚动轴承有限元模型

1.1 单元算法

显式算法以中心差分算法进行动态问题时域积分，其优点是无需考虑收敛性问题和解联立方程组。显式算法对于求解非线性问题、大变形问题和瞬态问题具有很好的效果。

本文所选Solid164单元采用单点积分算法，并用Lagrange列式进行求解计算，其列式的单元附着在材料上，可随材料的流动而产生单元网格变形，适合分析轴承这种变形较小的结构。

1.2 边界条件和材料参数

在实际工作中，轴承外圈一般嵌套在轴承座中不做运动。为模拟外圈在轴承座中的情况，将仿真模型中的轴承外圈外表面设置为刚性材料，并约束其所有自由度。保持架主要起引导滚动体运动，防止滚动体之间发生直接接触摩擦作用，外侧通常有端盖保护，为模拟端盖作用，将保持架X方向的平动自由度进行约束。内圈通常与刚性较大的轴连在一起，因此设置内圈内表面为刚性面，并约束内圈内表面X、Y、Z方向的平动自由度和Y、Z两个方向的转动自由度。根据试验轴承所在试验台的实际运转情况，设置轴承内圈内表面径向载荷为5kN，轴承内圈转速为150rad／s。轴承主体采用 GCr15钢，密度为7.85×10-6kg／mm3，弹性模量为206GPa，泊松比为0.3；保持架密度为7.85×10-6kg／mm3、弹性模量为196GPa、泊松比为0.24。

1.3 接触控制

轴承在运转过程中存在着3种接触，即滚动体与内圈滚道接触、滚动体与外圈滚道接触和滚动体与保持架兜空的接触。为提高接触的效率和接触的准确性，分析中采用指定主-从面的面-面接触方式，设置轴承内外圈滚道及保持架兜空为目标面，滚动体外表面为接触面，共设置24对接触。接触算法采用对称罚函数法［6］，其基本算法原理为：每一时刻检查从节点是否穿透主表面，若没有穿透则不作处理；如果存在穿透则在该从节点与被穿透主表面之间引入一个较大的界面接触力，即附加的法向接触力Fs：

式中，Δi为穿透量；Ev、Ai、Vi分别为主段所在单元的体积模量、主段面积和体积；f为接触刚度罚因子（缺省值为0.1）。

轴承运转过程中不可避免地存在着滚动摩擦和滑动摩擦，ANSYS／LS-DYNA的摩擦接触基于库仑公式，设摩擦因数为μ，静摩擦因数为μs，动摩擦因数为μd，则

式中，v为接触表面之间的相对速度；DC为衰减系数。

考虑摩擦的影响，将μ作适当放大，设置轴承在运转过程中滚动体与内外圈滚道的静摩擦因数为0.2，动摩擦因数为0.1；滚动体与保持架的静摩擦因数为0.1，动摩擦因数为0.05。

2 仿真结果分析

根据上述建立的模型及参数，将分析时间设置为0.12s，输出步数设为1000进行求解。

2.1 元件的应力分析

取轴承运转0.099 45s时刻的等效应力云图进行分析，如图2所示，从图中可以看出，轴承的最大应力发生在滚动体与内外圈接触的地方，且承载区应力要比非承载区应力大。保持架最大应力的发生位置是跟随兜空内滚动体最大应力的发生而发生的，但应力值较小。图3所示为各元件单元应力值的比较，从图中可明显看出，对于存在故障的外圈来说其应力值要远远大于内圈滚动体及保持架的应力值，外圈最大应力值约为450MPa，是内圈最大应力值的4倍，这与正常轴承滚动体最大应力值的情况存在明显不同。从图3还可看出各元件单元应力最大值的分布规律：外圈应力大于内圈应力，内圈应力大于滚动体应力、滚动体应力大于保持架应力，外圈最大等效应力值存在明显周期性，最大应力值的发生时刻要早于内圈最大应力值的发生时刻。

图2 0.099 45s时轴承轴向和保持架的等效应力云图

图3 轴承各元件单元等效应力

2.2 元件的振动响应

在轴承故障诊断中，通常是先利用故障检测系统将轴承的振动信号转换为数字信号，再利用信号分析技术进行分析，其中，最常用的是信号的FFT变换。现分别取外圈节点的X向、Y向、Z向（对应轴向、水平、竖直）位移、速度及加速度信号进行FFT变换分析。选择的节点位置如图1a所示，自故障位置起顺时针分别为节点1、节点2、节点3。图4a为接近故障位置的节点1的X向位移时频图，从中可明显看到周期性的振动冲击，但在频域图中并未找到与轴承外圈故障特征频率相近的频率成分，原因可能是X向为轴承旋转方向，在Z向负载的作用下，故障频率表现不明显。图4b、图4c分别为节点1的Y向和Z向位移时频图，从中可以看到明显的周期性振动冲击，结合图3分析可知，这种冲击是由轴承外圈周期性的应力集中造成的，因频域图中的72.93Hz与外圈故障特征频率73.1Hz极为接近，误差仅为0.2%，并存在明显的倍频成分，可以肯定为轴承外圈故障。图4b、图4c也存在着不同，图4b的时域信号仅存在负方向冲击，而图4c的时域信号为正负冲击交替出现，从幅值上看无论时域还是频域，图4b的信号幅值都比图4c的信号幅值大。以上分析表明，在检测位置的选择上，Y方向的信号比其他2个方向的信号所携带的信息要多，因此也更能体现故障特征。在节点速度响应分析中，能找到故障特征频率但是并不十分明显（图5a），而在图5b所示的加速度信号分析中未发现明显的故障特征频率，原因在于在数值的微分过程中存在着信息的丢失。

图4 节点1的X向、Y向和Z向位移时频图

图5 节点1的Y向振动响应时频图

所标识的节点2、节点3与节点1具有相似的情形，即Y方向的位移响应信号可明显地体现轴承的外圈故障，如图6所示。结合图4b及图6可看出，虽然三者均可明显地反映故障特征，但信号幅值存在明显差别，图6b（节点3）的幅值要高于图4a（节点1）、高于图6a（节点2）的幅值。Y方向的速度响应信号可找到特征频率成分但不明显，如图7a、图7b所示。从图7b（节点3）的速度信号中并没找到明显的故障特征信息，而从图5a（节点1）与图7a（节点2）中可找到故障特征频率值。在所取的图5b、图7c和图7d3个节点的加速度信号中均未找到轴承外圈故障的特征频率，因此，在分析故障信号时应尽可能地选择位移信号，位置的选择上应首先考虑Z向和Y向。

图6 节点2、3的Y向位移响应时频图