含局部故障的滚动轴承动力学建模及振动分析

2020-04-08李鸿飞

东北大学学报（自然科学版） 2020年3期

马辉，李鸿飞，俞昆，曾劲

(1.东北大学机械工程与自动化学院，辽宁沈阳 110819; 2.东北大学航空动力装备振动及控制教育部重点实验室，辽宁沈阳 110819)

滚动轴承是广泛应用于旋转机械中的重要支承零件，一旦发生故障将有可能影响整个设备的正常运行，造成重大的经济损失.在滚动轴承漫长的服役期内，有可能出现压痕、剥落、磨损、点蚀等各种局部故障，对含局部故障的滚动轴承进行动力学建模及振动分析具有重要的理论价值和工程意义，因而得到了国内外学者的广泛关注.

Zhang等[1]建立了滚动轴承显式有限元模型，提供了滚子滚过故障区的内部视角.李国超等[2]建立了滚动轴承的三维有限元模型，将局部缺陷引入健康轴承，研究了其运转及振动特点.Ding等[3]利用显式有限元方法研究了外圈局部故障尺寸对球轴承振动响应的影响.Liu等[4]基于滚动轴承的二维有限元模型研究了局部故障形状对球轴承振动波形的影响.Edwin[5]基于商用软件Algor建立了滚动轴承的二维模型，研究了不同载荷类型下滚子与滚道之间的接触特性.Singh等[6-7]建立了含外圈线剥落故障的二维滚动轴承显式有限元模型，通过多种信号处理方法详细阐述了低速重载下滚子滚过故障区的过程.Liu等[8]在此模型的基础上研究了不同缺陷边缘对故障轴承接触力和振动水平的影响.Liu等[9]对平面应力、平面应变、三维实体单元所建轴承模型进行了接触变形、接触应力方面的对比.Edwin等[10]基于Algor建立了滚动轴承的二维显式有限元模型，研究了滚子与滚道之间的打滑现象.张宇等[11]用显式有限元方法计算了球轴承在不同径向载荷下由加速到恒速过程的动力学响应.常斌全等[12]将滚动体经过内圈缺陷的过程进行了充分细化，建立了含局部缺陷的轴承的动力学模型.刘静等[13]考虑轴承滚子的凸度，研究了位移激励形式和局部故障尺寸对轴承振动特性的影响规律.

上述研究大致可分为两类：一类是考虑滚子与滚道之间的非线性接触而构建的集中质量模型，这类模型方程数量少，求解快，但是无法充分考虑轴承各零件的柔性，且只能得到轴承的振动响应；另一类是显式有限元模型，这类模型更多关注轴承外圈故障的情况，对轴承内圈故障的情况关注较少.有鉴于此，本文将首先建立一种圆柱滚子轴承的二维显式有限元动力学模型.然后在该模型中分别植入内外圈故障，分析故障区边缘单元的等效应力，滚子滚过内外圈局部缺陷的过程，将仿真所得振动响应与实验测试结果进行对比，以验证模型的有效性.最后研究故障区平滑程度对轴承振动特性的影响.

1 健康滚动轴承动力学模型

1.1 滚动轴承参数

为了方便实验轴承的加工，本文研究过程中选用两种型号的圆柱滚子轴承：1)对于健康和外圈故障的情况选用外圈无挡边的N205EM圆柱滚子轴承，具体参数见表1；2)对于内圈故障的情况选用内圈无挡边的NU205EM圆柱滚子轴承，具体参数见表2.

表1 N205EM轴承零部件几何尺寸

表2 NU205EM轴承零部件几何尺寸

若轴承外圈嵌在轴承座里固定不动，内圈随转轴同步转动，则恒转速下轴承的故障特征频率计算式为[14]

(1)

式中：fo为外圈故障特征频率；fi为内圈故障特征频率；Z为滚子数量；fs为转频；d为滚子直径；D为轴承节圆直径；α为轴承压力角，本文所用的圆柱滚子轴承的压力角为0°.

1.2 滚动轴承显式有限元建模

基于ANSYS/LS-DYNA建立N205EM轴承的二维显式有限元模型，如图1所示.整个轴承由Shell163单元建模，该单元由4个节点组成.轴承内圈、外圈、滚子的材料为轴承钢，密度为7 850 kg/m3，弹性模量为206 GPa，泊松比为0.3，保持架材料为黄铜，密度为8 920 kg/m3，弹性模量为100 GPa，泊松比为0.36.为避免初始穿透并且更加贴近真实的轴承，滚子与轴承外圈、内圈、保持架之间均留有一定间隙.轴承各零件之间的相互作用采用ASS2D算法模拟，阻尼比为2%.将轴承内圈最内侧的一层单元设为刚体以施加转速，约束轴承外圈最外层节点的所有自由度，以模拟轴承安装在轴承座中的情形.

图1 二维滚动轴承有限元模型

1.3 滚动轴承动力学模型的求解[15]

在LS-DYNA中，本文所建的有限元模型采用中心差分法求解.已知离散化的结构动力方程为

(2)

(3)

(4)

其中，c为声速，Ls为特征尺度：

(5)

(6)

式中：E为材料的弹性模量；ρ为材料密度；υ为泊松比；As是壳单元的面积;Li(i=1，2，3，4)是壳单元各边的长度.

1.4 滚动轴承有限元模型网格收敛性分析

在所建立的健康轴承有限元模型的内圈施加转速，前0.05 s为斜坡加载，0.05 s后为恒速700 r/min.保持架的理论转速由下式计算:

(7)

不同网格尺寸下，保持架转速与理论值的对比见图2.网格尺寸为0.5和0.25 mm时，保持架转速与理论值相符，从运动学的角度证明了本文模型的有效性.兼顾计算效率与求解精度，本文选取网格尺寸为0.5 mm，整个模型由10 363个单元，11 441个节点组成.

2 含局部故障的滚动轴承动力学建模

2.1 局部故障区域模拟

本文采用圆周矩形来模拟出现在轴承滚道上的通槽剥落.定义轴承外圈故障如图3a所示，其中α1，α2分别为故障区内外两圆弧所对应的圆心角，d1为故障深度，其局部有限元模型见图3b.定义轴承内圈故障如图3c所示，其中α3，α4分别为故障区内外两圆弧所对应的圆心角，d2为故障深度，其局部有限元模型如图3d所示.

图2 轴承保持架转速变化情况

图3 轴承故障示意图

2.2 滚子经过故障区域时的单元等效应力分析

对于外圈故障，滚子由图3b中的A单元滚入故障区，由B单元滚出故障区，在转速为700 r/min 时，提取A，B两单元的等效应力时程图，如图4a所示，包络分析分别见图4b和图4c.可知单元A，B所受等效应力以外圈故障特征频率的倒数为时间间隔产生周期性变化，并且总体上看B单元的应力水平比A单元高，说明滚子滚出故障区比滚入故障区会产生更加激烈的冲击和碰撞.

对于内圈故障的情况，滚子先经过图3d中的C单元，再经过D单元，在转速为2 000 r/min时，C，D两单元的等效应力时程图见图5a，其包络谱见图5b和图5c.谱图中的频率成分分别为转频fs及其倍频2fs、内圈故障特征频率fi及其倍频2fi，以及故障特征频率与转频的组合频率.通过计算相应时间段内有效值的方法来比较C，D两单元的应力水平，D单元略高，说明滚子滚出故障区会产生比滚入故障区产生更加激烈的冲击和碰撞，这一点与外圈故障时的情形一致.

图4 两单元等效应力(外圈故障)

图5 两单元等效应力(内圈故障)

2.3 滚子经过故障区域时的接触应力分析

图6所示为滚子滚过外圈故障区的应力云图.可分为三个过程：1) 滚子滚入故障区，在滚子与内外圈接触的地方应力非常大，随后滚子与故障区左侧边缘脱离，进入故障区；2) 滚子陷入故障区时与内外圈几乎不接触，作用在轴承内圈上的载荷全由其他滚子分担；3) 随着内圈的转动，陷在故障区的滚子被保持架带到故障区的右侧边缘，“楔入”内外圈之间的夹缝，与故障尖边发生碰撞，伴有比滚入故障区时更大的应力.

图6 滚子滚过外圈故障区

图7为滚子滚过内圈故障区的应力云图，具体的过程与外圈故障的情形相似，不再赘述.

2.4 含局部故障的滚动轴承动力学模型实验验证

本文用到的滚动轴承-转子系统实验台和待测轴承见图8.转轴左端由健康轴承支撑，右端由待测轴承支撑.实验过程中，采样频率设为20 kHz.按图3b和图3d分别在N205EM轴承外圈和NU205EM轴承内圈上加工圆周矩形的凹槽.

图7 滚子滚过内圈故障区

图8 实验台、故障轴承

对于外圈故障，当转速为700 r/min时，距轴承外圈外圆1.5 mm，周向位置0°，90°，180°，270°(以x轴正方向为0°，逆时针为正)的4个节点的加速度信号如图9a～图9d所示.0°，180°方向上的节点加速度信号比较相似，而位于载荷区的270°的节点冲击比较明显.因此本文选取载荷区节点的加速度与实验进行对比，如图10所示.实验与仿真的时域加速度振动信号都很好地反映了由于外圈局部故障的存在而导致的冲击现象，且本文所用有限元模型充分考虑了轴承各零件的柔性和惯性，所得时域信号的波形与实验吻合较好.加速度信号的包络谱中轴承外圈故障特征频率及其倍频清晰可见.仿真、实验、理论(根据式(1)计算)所得的轴承外圈故障特征频率对比见表3，由此说明本文所建的有限元模型能够较准确地模拟外圈故障时的滚动轴承振动特性.

表3 不同方法所得外圈故障特征频率

图9 周向节点加速度对比

图10 仿真、实验对比图(外圈故障)

对于内圈故障，当转速为2 000 r/min时，仿真结果与实验实测结果的对比见图11.轴承内圈存在局部故障时，由于受到转频的调制作用，包络谱中可以看到由于故障导致的内圈故障特征频率fi及其倍频2fi，以及由于转频调制而出现的故障特征频率与转频fs的组合频率.仿真、实验、理论(根据式(1)计算)所得的轴承内圈故障特征频率对比见表4，二者相近，由此证明了本文所建的有限元模型也可以较为准确地模拟内圈故障时的滚动轴承振动特性.