2020, Vol. 41
2. 东北大学 航空动力装备振动及控制教育部重点实验室，辽宁 沈阳 110819
2. Key Laboratory of Vibration and Control of Aero-Propulsion System, Ministry of Education, Northeastern University, Shenyang 110819, China
滚动轴承是广泛应用于旋转机械中的重要支承零件，一旦发生故障将有可能影响整个设备的正常运行，造成重大的经济损失.在滚动轴承漫长的服役期内，有可能出现压痕、剥落、磨损、点蚀等各种局部故障，对含局部故障的滚动轴承进行动力学建模及振动分析具有重要的理论价值和工程意义，因而得到了国内外学者的广泛关注.
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轴承零部件几何尺寸 Table 1 Dimensions of N205EM components |
|
|
表 2 NU205EM轴承零部件几何尺寸 Table 2 Dimensions of NU205EM components |
若轴承外圈嵌在轴承座里固定不动，内圈随转轴同步转动，则恒转速下轴承的故障特征频率计算式为[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 二维滚动轴承有限元模型 Fig.1 A meshed 2-D FE model of rolling bearing |
在LS-DYNA中，本文所建的有限元模型采用中心差分法求解.已知离散化的结构动力方程为
|
(2) |
式中：M为系统质量矩阵；C为阻尼矩阵；ẍ(t)，ẋ(t)分别为节点的加速度向量和速度向量；P(t)，F(t)，H(t)分别为荷载向量、内力向量和沙漏阻力向量.上述系统动力方程采用中心差分法求解的基本递推格式为
|
(3) |
式中：t(n-1)/2=(tn+tn-1)/2，t(n+1)/2=(tn+tn+1)/2，Δtn-1=tn-tn-1，Δtn=tn+1-tn；ẍ(t)，ẋ(t)，x(t)分别为对应时刻的节点加速度、速度和位置坐标向量.这是一种显式求解方法，其中每一时刻的积分步长由稳定性条件控制，对于本文中的壳单元，其表达式为
|
(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 轴承保持架转速变化情况 Fig.2 Variation of angular velocity of bearing cage |
本文采用圆周矩形来模拟出现在轴承滚道上的通槽剥落.定义轴承外圈故障如图 3a所示，其中α1，α2分别为故障区内外两圆弧所对应的圆心角，d1为故障深度，其局部有限元模型见图 3b.定义轴承内圈故障如图 3c所示，其中α3，α4分别为故障区内外两圆弧所对应的圆心角，d2为故障深度，其局部有限元模型如图 3d所示.
|
图 3 轴承故障示意图 Fig.3 Bearing fault diagram (a)—外圈故障示意图；(b)—外圈故障有限元图；(c)—内圈故障示意图；(d)—内圈故障有限元图. |
对于外圈故障，滚子由图 3b中的A单元滚入故障区，由B单元滚出故障区，在转速为700 r/min时，提取A，B两单元的等效应力时程图，如图 4a所示，包络分析分别见图 4b和图 4c.可知单元A，B所受等效应力以外圈故障特征频率的倒数为时间间隔产生周期性变化，并且总体上看B单元的应力水平比A单元高，说明滚子滚出故障区比滚入故障区会产生更加激烈的冲击和碰撞.
|
图 4 两单元等效应力(外圈故障) Fig.4 Equivalent stress of two elements (outer ring fault) (a)—A, B两单元等效应力时程图；(b)—A单元等效应力包络谱；(c)—B单元等效应力包络谱. |
对于内圈故障的情况，滚子先经过图 3d中的C单元，再经过D单元，在转速为2 000 r/min时，C，D两单元的等效应力时程图见图 5a，其包络谱见图 5b和图 5c.谱图中的频率成分分别为转频fs及其倍频2 fs、内圈故障特征频率fi及其倍频2fi，以及故障特征频率与转频的组合频率.通过计算相应时间段内有效值的方法来比较C，D两单元的应力水平，D单元略高，说明滚子滚出故障区会产生比滚入故障区产生更加激烈的冲击和碰撞，这一点与外圈故障时的情形一致.
|
图 5 两单元等效应力(内圈故障) Fig.5 Equivalent stress of two elements (inner race fault) (a)—C, D单元等效应力时程图；(b)—C单元等效应力包络谱；(c)—D单元等效应力包络谱. |
图 6所示为滚子滚过外圈故障区的应力云图.可分为三个过程：1)滚子滚入故障区，在滚子与内外圈接触的地方应力非常大，随后滚子与故障区左侧边缘脱离，进入故障区；2)滚子陷入故障区时与内外圈几乎不接触，作用在轴承内圈上的载荷全由其他滚子分担；3)随着内圈的转动，陷在故障区的滚子被保持架带到故障区的右侧边缘，“楔入”内外圈之间的夹缝，与故障尖边发生碰撞，伴有比滚入故障区时更大的应力.
|
图 6 滚子滚过外圈故障区 Fig.6 Roller rolls over the outerring fault zone (a)—滚入故障区；(b)—陷落故障区；(c)—滚出故障区. |
图 7为滚子滚过内圈故障区的应力云图，具体的过程与外圈故障的情形相似，不再赘述.
|
图 7 滚子滚过内圈故障区 Fig.7 Roller rolls over the innerring fault zone (a)—脱离；(b)—陷落；(c)—相撞. |
本文用到的滚动轴承-转子系统实验台和待测轴承见图 8.转轴左端由健康轴承支撑，右端由待测轴承支撑.实验过程中，采样频率设为20 kHz.按图 3b和图 3d分别在N205EM轴承外圈和NU205EM轴承内圈上加工圆周矩形的凹槽.
|
图 8 实验台、故障轴承 Fig.8 Test rig and the fault bearings (a)—滚动轴承-转子系统实验台；(b)—外圈故障轴承；(c)—内圈故障轴承. |
对于外圈故障，当转速为700 r/min时，距轴承外圈外圆1.5 mm，周向位置0°，90°，180°，270°(以x轴正方向为0°，逆时针为正)的4个节点的加速度信号如图 9a~图 9d所示.0°，180°方向上的节点加速度信号比较相似，而位于载荷区的270°的节点冲击比较明显.因此本文选取载荷区节点的加速度与实验进行对比，如图 10所示.实验与仿真的时域加速度振动信号都很好地反映了由于外圈局部故障的存在而导致的冲击现象，且本文所用有限元模型充分考虑了轴承各零件的柔性和惯性，所得时域信号的波形与实验吻合较好.加速度信号的包络谱中轴承外圈故障特征频率及其倍频清晰可见.仿真、实验、理论(根据式(1)计算)所得的轴承外圈故障特征频率对比见表 3，由此说明本文所建的有限元模型能够较准确地模拟外圈故障时的滚动轴承振动特性.
|
图 9 周向节点加速度对比 Fig.9 Nodes acceleration comparison along circumferential direction (a)—0°；(b)—90°；(c)—180°；(d)—270°. |
|
图 10 仿真、实验对比图(外圈故障) Fig.10 Comparison between simulation and experiment (outer race fault) (a)—仿真加速度时域图；(b)—仿真加速度包络谱；(c)—实验加速度时域图；(d)—实验加速度包络谱. |
|
|
表 3 不同方法所得外圈故障特征频率 Table 3 Comparison of outer race fault frequency using different methods |
对于内圈故障，当转速为2 000 r/min时，仿真结果与实验实测结果的对比见图 11.轴承内圈存在局部故障时，由于受到转频的调制作用，包络谱中可以看到由于故障导致的内圈故障特征频率fi及其倍频2fi，以及由于转频调制而出现的故障特征频率与转频fs的组合频率.仿真、实验、理论(根据式(1)计算)所得的轴承内圈故障特征频率对比见表 4，二者相近，由此证明了本文所建的有限元模型也可以较为准确地模拟内圈故障时的滚动轴承振动特性.
|
图 11 仿真、实验对比图(内圈故障) Fig.11 Comparison between simulation and experiment (inner race fault) (a)—仿真加速度时域图；(b)—仿真加速度包络谱；(c)—实验加速度时域图；(d)—实验加速度包络谱. |
|
|
表 4 不同方法所得内圈故障特征频率 Table 4 Comparison of inner race fault frequency using different methods |
对于外圈故障的情况，假设故障深度d1=0.5 mm，圆心角α1=5°，随α2增大，故障轮廓越来越平滑.提取不同α2下的轴承振动加速度响应并对其进行包络分析，如图 12所示.可见，故障轮廓越平滑，滚子滚出故障区楔进内外圈的夹缝越容易，冲击越小.
|
图 12 不同角度下的节点加速度(外圈故障) Fig.12 Nodes acceleration under different angles(outer ring fault) (a)—不同圆心角下的加速度；(b)—不同圆心角下的包络谱. |
对于内圈故障的情况，假设故障深度d2=0.5 mm，圆心角α3=5°，提取不同α4下的轴承振动加速度响应并对其进行包络分析，如图 13所示，幅值总体上呈现先增大后减小的趋势.
|
图 13 不同角度下的节点加速度(内圈故障) Fig.13 Nodes acceleration under differrnt angles(inner ring fault) (a)—不同圆心角下的加速度；(b)—不同圆心角下的包络谱. |
1) 无论外圈故障还是内圈故障，滚子滚出故障区域时，与故障区域边缘接触，产生更强烈的冲击和碰撞；
2)对外圈故障，随圆心角α2增大，轴承的故障轮廓越来越平滑，由此产生的振动响应越来越小；
3)内圈故障时，随圆心角α4增大，故障轮廓的变化先以宽度增加占主导，在圆心角超过某一临界值后，故障轮廓逐渐平滑，因此幅值先大后小.
| [1] |
Zhang Z N, Ding W M, Ma H F. Local stress analysis of a defective rolling bearing using an explicit dynamic method[J]. Advances in Mechanical Engineering, 2016, 8(12): 1-9. |
| [2] |
李国超, 彭炜, 李勇才, 等. 滚动轴承外圈故障的显式有限元动态仿真分析[J]. 中国机械工程, 2012, 23(23): 2825-2829. (Li Guo-chao, Peng Wei, Li Yong-cai, et al. Simulation and dynamic analysis of outer ring fault on rolling bearings using explicit finite element method[J]. China Mechanical Engineering, 2012, 23(23): 2825-2829. DOI:10.3969/j.issn.1004-132X.2012.23.011) |
| [3] |
Ding W M, Zhang Z N, Zhao F G.Vibration response of ball bearings with different defect sizes in the outer raceway: simulation with a 3-D finite element model[C]//The 14th IFToMM World Congress.Taipei, 2015: 534-537.
|
| [4] |
Liu J, Shao Y M, Zuo M J. The effects of the shape of localized defect in ball bearings on the vibration waveform[J]. Proceedings of the Institution of Mechanical Engineers Part K:Journal of Multi-body Dynamics, 2013, 227(3): 261-274. |
| [5] |
Edwin L J. Numerical model to study of contact force in a cylindrical roller bearing with technical mechanical event simulation[J]. Journal of Mechanical Engineering & Automation, 2011, 1(1): 1-7. |
| [6] |
Singh S, Köpke U G, Howard C Q, et al. Analyses of contact forces and vibration response for a defective rolling element bearing using an explicit dynamics finite element model[J]. Journal of Sound & Vibration, 2014, 333(21): 5356-5377. |
| [7] |
Singh S, Howard C Q, Hansen C H, et al. Analytical validation of an explicit finite element model of a rolling element bearing with a localized line spall[J]. Journal of Sound & Vibration, 2018, 416: 94-110. |
| [8] |
Liu J, Shao Y M. A numerical investigation of effects of defect edge discontinuities on contact forces and vibrations for a defective roller bearing[J]. Proceedings of the Institution of Mechanical Engineers Part K:Journal of Multi-body Dynamics, 2015, 230(4): 387-400. |
| [9] |
Liu J, Wu H, Shao Y M. The influence of the raceway thickness on the dynamic performances of a roller bearing[J]. Journal of Strain Analysis for Engineering Design, 2017, 52(2): 1-9. |
| [10] |
Edwin L J, Jesus M A, Vicente D L. A study of sliding between rollers and races in a roller bearing with a numerical model for mechanical event simulations[J]. Tribology International, 2010, 43(11): 2175-2182. |
| [11] |
张宇, 谢里阳, 吴宁祥, 等. 深沟球轴承加速-恒速过程动态特性有限元分析[J]. 东北大学学报(自然科学版), 2012, 33(1): 103-107. (Zhang Yu, Xie Li-yang, Wu Ning-xiang, et al. Finite element analysis for dynamic characteristic of deep-groove ball bearing operating from acceleration to constant speed[J]. Journal of Northeastern University(Natural Science), 2012, 33(1): 103-107.) |
| [12] |
常斌全, 剡昌锋, 苑浩, 等. 多事件激励的滚动轴承动力学建模[J]. 振动与冲击, 2018, 37(17): 16-24. (Chang Bin-quan, Yan Chang-feng, Yuan Hao, et al. Dynamic modeling for rolling bearings under multi-event excitation[J]. Journal of Sound & Vibration, 2018, 37(17): 16-24.) |
| [13] |
刘静, 邵毅敏, 秦晓猛, 等. 基于非理想Hertz线接触特性的圆柱滚子轴承局部故障动力学建模[J]. 机械工程学报, 2014, 50(1): 91-97. (Liu Jing, Shao Yi-min, Qin Xiao-meng, et al. Dynamic modeling on localized defect of cylindrical roller bearing based on non-Hertz line contact characteristics[J]. Journal of Mechanical Engineering, 2014, 50(1): 91-97.) |
| [14] |
Harris T, Kotzalas M.Rolling bearing analysis-essential concepts of bearing technology[M].5th ed.New York: Taylor and Francis, 2007: 220-250.
|
| [15] |
Hallquist J O. LS-DYNA theory manual[M]. California: Livermore Software Technology Corporation, 2006: 137-143.
|