考虑滚动轴承故障处塑性变形的有限元建模与动力学特性分析

涂文兵,杨本梦,杨锦雯,梁杰,梁才雪

(华东交通大学机电与车辆工程学院,330013,南昌)

随着机械装备向着高速、重载、高精度方向发展,对滚动轴承的动态服役性能的要求越来越高。目前,轴承在运行过程中出现故障的现象时有发生,这对机械装备高效运行产生了很严重的影响,因此对滚动轴承故障的动力学特性进行研究具有重要理论意义和工程应用价值。

滚动轴承局部故障的动力学建模是研究轴承故障机理和动力学特性的常用方法。Patil等[1]基于赫兹接触变形理论,将滚道和滚动体之间的接触力考虑成非线性接触弹簧,将局部故障激励考虑为半正弦位移激励函数,建立了滚动轴承的故障动力学模型,在此模型基础上研究了缺陷尺寸和位置对轴承振动特性的影响。东亚斌等[2]考虑了缺陷尺寸、载荷区和滚动体上点蚀缺陷滚过内外圈时的弹性接触力变化,建立了单一局部故障的滚动轴承模型,分析了当滚动体具有单一局部缺陷时轴承振动信号的成分。

上述研究大多数是采用恒值位移激励或恒定力激励代替局部故障产生的冲击激励,然而由于局部故障激励与轴承故障尺寸、轴承转速等条件密切相关,单一简化激励函数难以准确反映滚动体与故障的实际接触情况。

Liu等[3]基于分段函数和赫兹接触机制的耦合,建立了与球轴承滚道表面局部尺寸相关的球轴承故障动力学模型,提出了一种动态模拟方法,可以为不同缺陷尺寸的滚动体与滚道之间的接触提供更接近真实的脉冲。Shao等[4]提供了一种模拟局部表面缺陷对圆柱滚子轴承振动响应影响的新方法,建立了一个考虑时变的挠度激励和时变的接触刚度激励的二自由度(DOF)的动态模型。常斌全等[5]考虑径向载荷作用下轴承缺陷引起的时变位移对动力学响应特性的影响,将滚动体过缺陷的接触过程细化为进入、冲击、离开和载荷补偿4个过程,建立了内圈表面局部缺陷的滚动轴承二自由度动力学模型,研究了缺陷尺寸和转速对轴承系统振动加速度幅值的影响规律。罗茂林等[6]基于赫兹接触理论,以外圈滚道表面存在剥落的球轴承为研究对象,建立了考虑冲击力的球轴承外圈剥落双冲击现象动力学模型。Yang等[7]利用脉冲序列模型和转子-轴承-机壳系统中的非线性多体动态轴承模型对轴承故障特征的原因和影响因素进行了解释,分析了缺陷位置和随机波动的影响。李昊泽等[8]基于刚性套圈假设和赫兹接触理论,考虑时变位移激励函数,建立了滚动轴承外圈局部故障动力学模型,研究了故障尺寸对一维及多维轴承振动特征的影响规律。胡爱军等[9]采用五自由度动力学模型,分析了滚动轴承外圈故障数、故障间隔和载荷分布对故障特征的影响。Patra等[10]考虑到内外圈滚道和滚子等主要部件上特定尺寸的剥落缺陷,进行了定性分析。张慧玲等[11]考虑油膜和滚动体滑动等非线性因素,建立了包含局部故障的滚动轴承系统五自由度动力学模型。李志农等[12]建立了含滚动体点蚀故障的滚动轴承动力学渐变模型,并给出了渐变释放的变形量的计算方法。王震等[13]提出了非规则轴承故障的表征方法,建立了轴承内外圈非规则故障模型。揭示了故障的周向宽度和深度对系统振动的影响规律。王凯等[14]针对径向载荷作用下复合故障激励的滚动轴承振动机理的复杂性问题,提出了四自由度复合故障深沟球轴承动力学模型,分析了滚动轴承复合故障对系统动力学振动响应的影响。Zhang等[15]建立了基于时变位移的滚动轴承四自由度动态模型,研究了不同条件下滚动轴承局部缺陷的振动特性。通过时域和频域的实验振动信号对建立的模型进行了验证。

上述模型虽然获取了轴承的动态响应,但是大多数模型都是将轴承的实际故障用位移激励来代替,且对轴承零部件作了刚性假设,与实际情况存在一定的偏差。轴承在实际运行过程中会产生一定的变形,这在传统的动力学建模中难以考虑。随着轴承研究的进一步推进,以及计算机技术的不断发展,有限元仿真分析(FEA)为分析模拟轴承实际运行过程中的动力学特性问题提供了有效的帮助。林腾蛟等[16]在ANSYS/LSDYNA环境中建立了深沟球轴承的三维有限元模型并进行计算,将有限元计算的结果和赫兹理论计算的结果进行比较,两种方法所得结果基本一致。

考虑轴与轴承座之间的相互影响,刘静等[17-18]考虑轴承各元件之间的接触与摩擦的影响以及剥落故障对滚子与故障周边区域之间的接触特性和轴承振动响应特征的影响,分别建立含局部故障的轴承-轴承座系统有限元动力学模型和考虑故障过渡区的圆柱滚子轴承动力学模型,分析不同局部故障尺寸变化对其振动特征的影响规律。马辉等[19]建立了正常和含局部剥落故障的圆柱滚子轴承的二维有限元模型,分析了故障区边缘单元的等效应力和滚子滚过内外圈局部缺陷的过程,研究了局部故障区域平滑程度对滚动轴承振动特性的影响,结果表明,故障区域越平滑,振动响应越小。Tu等[20]建立了外滚道三类故障的动态非线性有限元模型,将局部故障表征为滚道表面的矩形凹槽,分析了滚动体动态接触力和运动特性,揭示了故障相关力随滚动体轴承故障演变而产生振动响应的机制。杨子臻等[21]运用LS-DYNA建立转子-轴承系统二维显式动力学有限元模型,分析了故障耦合状态下轴承振动加速度、转子轴心运动轨迹,探明了转子不平衡-轴承缺陷耦合故障冲击响应特征及转子轴心运动轨迹突变机理。

有限元仿真分析方法可以很好地考虑轴承各部件的变形,更加符合轴承的实际运行状况,但是在以往有限元仿真分析中将轴承的变形理想化为弹性变形,然而轴承在实际运行过程中,尤其是故障产生后的局部接触容易引发应力集中现象,导致局部出现明显的塑性变形,对滚动的接触状态和运动特性产生重要影响,进而影响故障轴承的动力学特性。虽然Howell等[22]运用了一种纯滚动的线接触的二维有限元模型,分析了持续循环变形、蚀刻带和裂纹之间的联系。Liu等[23]在有限元模型中采用了双线性随动强化模型,建立了边缘塑性变形与剥落边缘轮廓之间的关系,但是都没有考虑塑性变形对振动的影响。

本文以NU306轴承为研究对象,在ANSYS LS-DYNA环境下建立了考虑滚动轴承故障塑性的3种动力学模型,分别为弹/弹模型、塑/塑模型、弹/塑模型,选用塑性随动强化模型和线弹性材料模型分别计算材料的塑性变形和弹性变形,通过将3种动力学模型仿真结果与实验结果进行对比,揭示了故障状态下轴承产生塑性变形与响应信号的内在联系,为实际运行中塑性变形对故障轴承动力学特性的影响分析提供理论依据。

1 有限元模型的建立

1.1 轴承的主要结构参数

本文以NU306圆柱滚子轴承作为研究对象,建立了轴承-轴承座系统的显式动力学有限元模型,其主要参数如表1所示。

表1 NU306轴承主要参数

1.2 单元类型及网格划分

圆柱滚子轴承只能承受径向载荷,且结构具有一定的轴向对称性,考虑到内部运动关系复杂和非线性接触等因素,在满足计算精度的情况下尽量节约时间,因此本文采用二维平面单元来建立轴承-轴承座系统的有限元模型。由于四边形网格计算时更加容易收敛,计算结果更加稳定,本文采用Shell 181二维应变单元类型对模型进行四边形网格划分。网格划分必须保证滚动体表面和滚道在模拟过程中始终保持接触,因此需要设置一定的计算精度,为了满足以上要求,本文中模型设置了大小分别为1 mm、0.75 mm、0.5 mm、0.25 mm和0.125 mm的网格尺寸进行试运算。对于采用1 mm和0.75 mm网格元素的有限元模型,发现在模拟过程中,滚动体和滚道之间不能保持连续接触,影响滚动体上载荷分布;对于采用0.25 mm和0.125 mm的尺寸进行网格划分的有限元模型,上述条件得到了满足,从而能够对轴承运动学特性进行准确获取,由于上述两种网格仿真所得结果非常接近,因此,综合考虑计算机硬件和计算效率,为了减少CPU的运行时间而不影响数值解的准确性,选择0.25 mm的网格尺寸模型。为了更好地模拟两元件之间的接触,在滚动体与内外圈及保持架接触部位进行了网格细化,最终划分单元数为85 028个,建模时采用宽度和深度都为1 mm的贯穿式矩形坑来模拟故障,带故障的滚动轴承有限元模型如图1所示。

图1 带故障的滚动轴承有限元模型Fig.1 Finite element model of rolling bearings with faults

1.3 材料模型及参数

考虑到产生塑性变形的不同位置,为全面分析塑性变形的影响,本文建立了故障轴承的3种有限元模型,分别为弹/弹模型(滚动体和外圈采用线弹性材料本构模型)、弹/塑模型(滚动体采用线弹性材料本构模型,外圈采用塑性随动强化模型)、塑/塑模型(将滚动体和外圈都采用塑性随动强化模型)。保持架采用黄铜,圆柱滚子和内外圈及轴承座材料均采用GCr15钢,其材料参数如表2所示,塑性随动强化模型参数方程如下式[24-25]

(1)

(2)

在定义该模型的时候需要输入参数包括弹性模量(EX)、密度(DENS)、泊松比(NUXY)、屈服强度和切变模量,这些参数[26]在表2中给出。

表2 NU306滚动轴承及轴承座模型材料参数

1.4 接触类型及边界条件

LS-DYNA中通常采用罚函数进行接触结算,通过增大接触力罚常数来减少穿透量,避免初始穿透影响计算结果的有效性。另外,摩擦力是轴承运行的重要影响因素,LS-DYNA中基于库伦公式计算摩擦力,其摩擦因数μ计算方式如下式

μ=μd+(μs-μd)e-c|v|

(3)

式中:v表示接触面间的相对速度;c为衰减系数;为了模拟轴承在油润滑方式下的运动情况,设置静摩擦因数μs=0.15,动摩擦因数μd=0.002[20]。相对点面接触来说,面面接触更能对高、低阶单元的大滑动和大变形刚度进行更好的协调,因此选用二维接触中的自动面面接触,定义12个滚动体为接触面,内、外圈滚道及保持架为目标面,建立36组接触对来模拟滚动体与各元件间的高度非线性接触行为。

在动态响应计算中,经常会出现沙漏现象,导致模型显示不符合实际的非常奇怪的网格,因此在网格划分时尽量均匀细化,避免单点载荷,尽量使用全积分运算,以及通过控制线型和二次系数来增大模型的体积黏性等,达到控制沙漏的目的。综合上述方法,设置沙漏系数为0.1[26],该值超过0.15将会引起计算不稳定,体积黏度系数使用默认值1.5,线性体积黏度系数使用默认值0.06。

对于边界条件的设置,在实际工作中,滚动轴承外圈一般全约束在轴承座内,轴承内圈绕Z轴旋转,内圈转速为2 100 r/min。此外NU306轴承正常运行时会承受一定径向力的作用,因此本文在轴承内圈内表面对Y方向的负方向施加3 000 N的径向力;为了避免轴承运转同时突然施加载荷对正常运转轴承的动态响应的影响,先施加时长0.02 s的径向载荷,待载荷稳定后再施加转速[27]。

2 振动特性分析

为了验证上述故障滚动轴承有限元模型的有效性,本文搭建了如图2所示的滚动轴承实验台。为保证实验运行条件与仿真模型运行条件相一致,在测试轴承内圈施加3 000 N的径向力,并通过伺服电机调整轴承内圈转速为2 100 r/min;为了保证测试轴承与仿真模型的一致性,采用电火花法在测试轴承外圈内表面加工宽度和深度都为1 mm的矩形故障。

图2 滚动轴承实验台Fig.2 Rolling bearing test bench and outer ring failure

提取滚动轴承外圈故障附近63 696节点的仿真与实验加速度时域信号如图3所示,特征频率及加速度幅值如表3所示。

图3 仿真与实验加速度时域信号Fig.3 Simulation and experimental acceleration time domain signals

当滚动体经过轴承外圈故障区域时会引起强烈振动,这种振动在加速度时域信号中出现周期性的冲击成分。对比图3中时域波形可以看出,仿真所得加速度时域波形与实验所得结果具有高度的一致性,实验加速度信号中相邻冲击信号周期为6.08 ms,仿真加速度信号中相邻冲击信号周期为6.10 ms,相对实验所测得结果误差仅为0.328%,有误差的原因可能是实际运行过程中滚动体与保持架之间的不稳定接触造成的。对比表3中仿真计算所得滚动体通过外圈故障的通过频率166.5 Hz,同实验所得轴承外圈故障特征频率164.7 Hz比较,特征频率基本一致,通过下式计算滚动轴承外圈故障特征频率

表3 特征频率及加速度幅值

(4)

式中:n为内圈转速;α为接触角。将模型参数代入得到滚动体通过外圈故障的通过频率165.15 Hz,同外圈故障特征频率的仿真结果166.5 Hz及实验结果164.7 Hz比较,得出其值具有高度的一致性,进一步说明该有限元模型对用于故障轴承的动力学分析是可行的。

对比表3中仿真与实验信号幅值可知,弹/弹模型加速度幅值远高于实验所测得的加速度幅值,而弹/塑模型及塑/塑模型所得加速度幅值与实验所得结果更加接近。根据表4中时域指标数据可知,弹/弹模型的振动加速度均方根1.557 6×105mm·s-2、峰峰值1.57×106mm·s-2约为弹/塑模型、塑/塑模型以及实验所得值的两倍多,而弹/塑模型、塑/塑模型与实际实验所测得结果相差不大。可能的原因是轴承在实际运行过程中,故障产生后局部接触容易引发应力集中,导致故障处产生塑性变形,滚动体与故障外圈接触力减小、振动信号减弱,其中弹/塑模型及塑/塑模型所得仿真结果与实验所得结果接近但略有差异的原因是由于实际实验时轴承润滑脂能够衰减一定的振动,导致实际测得振动信号减弱,另外仿真模型网格划分时会有多边形效应而引起的振动,导致仿真信号幅值以及时域指标数据略大于实验所得结果,由此再一次证明了本文所建立的轴承故障塑性动力学模型的有效性。

表4 故障信号的时域指标值

为了得到一个更加精确有效的故障塑性动力学模型,本文进一步分析了弹/塑模型和塑/塑模型之间的差异及对轴承振动特性的影响。由图3可知,考虑滚动体产生塑性变形的塑/塑模型与不考虑滚动体产生塑性变形的弹/塑模型所得加速度时域信号几乎一致,且从表3及表4中所得数据对比可知,弹/塑模型和塑/塑模型所得结果几乎一致,可能是因为在实际运行过程中,滚动体几乎不会产生塑性变形,因此在今后的建模过程中可以仅考虑故障处的塑性变形。

为了验证上述可能原因中滚动体及故障外圈塑性变形对轴承振动特性产生的影响,本文取图1所示模型中滚动体23 101单元以及故障处边缘61 332单元的应变进行分析,得出应变曲线如图4所示。

(a)滚动体上23 101单元

(b)故障处边缘61 332单元图4 滚动体及故障边缘应变曲线Fig.4 Strain curves of roller and fault edges

由图4(a)可以看出,滚动体变形量非常小,表现为塑/塑模型中,滚动体过故障时的应变值远小于滚动体为弹性材料时的应变值,当滚动体为塑性材料时其应变值会有塑性叠加现象,表现为进入承载区后,滚动体每次过故障时其应变值会有极小的突增并维持一定值,即使在非承载区,滚动体的应变也会维持在增加后的值不变,这是由于其变形不可恢复所导致,然而这种变形量都是极小的,几乎可以忽略不计;从图4(b)可知,故障产生后的局部接触容易引发应力集中现象,导致局部出现明显的塑性变形,因此塑/塑模型及弹/塑模型的应变值远大于弹/弹模型的应变值;除此之外,故障边缘61 332单元的应变值虽然略大于塑/塑模型该单元的应变值,但是两者数值相差极小,该结果说明滚动体塑性变形对轴承动态特性影响不大,实际运行过程中滚动体几乎不会产生塑性变形,为上述分析中仅考虑轴承故障处塑性动力学建模提供了有效的依据。

3 接触特性分析

为进一步揭示塑性变形所导致的故障轴承振动特性变化的本质原因,本节对滚动体过故障时接触特性进行分析。

3.1 滚动体与故障外圈的接触力分析

由于轴承中每一个滚动体与故障外圈接触力变化规律基本相同,因此本文以图1中7#滚动体与故障外圈的接触力曲线的单个周期进行分析,得出接触力曲线如图5所示。

图5 滚动体过故障时与外圈接触力变化曲线Fig.5 The contact force change curve of the roller over the out ring failure

分析图5及表5中3种模型的接触力曲线及接触力可知,该结果与文献[20]中所得出的结果具有一致性,当滚动体进入承载区时,其与故障外圈之间的接触力逐渐增大;刚进入故障时,会出现去接触力阶段(图5中ab、ac阶段),在此过程中,滚动体和外滚道之间的接触力由1 290 N开始逐渐减少,在图5中表现为7#滚动体通过故障时,其接触力迅速减小,对于弹/弹模型而言,其变形可恢复,故滚动体落入故障的深度不大,因此还有一定的接触力,其接触力减小至636 N;与弹/弹模型相比,弹/塑模型、塑/塑模型滚动体过故障时接触力减小的更多(图5中ac阶段),其值由1 290 N直接减小至0,这是因为含塑性材料模型故障变形量增大,导致滚动体与外圈完全失去接触。进一步分析可知弹/塑模型和塑/塑模型滚动体过故障时接触力下降程度几乎相同,根据上述应变分析可知,这是因为滚动体几乎不会产生塑性变形,因此建模时可以不用考虑滚动体塑性变形对轴承外圈故障塑性变形的影响。当滚动体准备从故障中离开时,弹/弹模型滚动体所受载荷逐渐增大(图5中bd阶段),而弹/塑模型、塑/塑模型滚动体与故障外圈接触力恢复时间略大于弹/弹模型所需时间(图5中cd阶段),在图5中表现为abd阶段包含在acd阶段内,这是因为弹/塑模型、塑/塑模型滚动体与外圈完全失去接触后,从故障中离开时需要重新与故障边缘接触,滚动体重新受力,滚动体与故障外圈接触力才能够逐渐增大。

表5 滚动体过故障时与外圈接触力

从图5中还可以看出,滚动体过故障时,承载区滚动体上的载荷会重新分配。当7#滚动体通过故障时,接触力减小,该滚动体上无法承受载荷,而相邻的6#以及8#滚动体接触力增大,这是由于滚道上的故障导致7#滚动体与外圈接触间隙增大,接触力减小,总载荷不变情况下,减小的这一部分载荷由其余滚动体承担,因此相邻两滚动体的接触力增大,进行载荷补偿;通过对表5中数值进行分析可知,含塑性材料模型7#滚动体相邻两滚动体接触力为976 N及以上,相对弹/弹模型过故障时相邻滚动体接触力939 N增大的更多,这是由于故障处塑性变形导致滚动体与滚道之间的间隙更大,接触力更小,相邻滚动体所需要承受的补偿量更大,图5中放大处清晰可见含塑性材料模型7#滚动体相邻两滚动体接触力增大的更多。

根据上述分析可知,在轴承实际运行过程中故障塑性变形会导致滚动体与故障处接触力显著减小,同时导致轴承中各滚动体受载荷不均,使得各滚动体运行状况不一致,导致轴承运行不稳定现象。

3.2 滚动体与保持架兜孔的接触力分析

滚动体与保持架兜孔之间接触力的变化会直接或者间接的对滚动体与滚道之间的接触力以及应变产生影响。由此,进一步探究3种模型中滚动体过故障时与保持架兜孔之间接触力的变化。滚动体通过故障期间与保持架兜孔的接触力曲线如图6所示。

图6 滚动体通过故障期间与保持架兜孔的接触力曲线Fig.6 The contact force curve between the roller and the cage pocket when the roller over the failure

理论上,承载区应该是滚动体推动保持架向前运动,接触力应全部为正值,但是对于有故障轴承,根据上述滚动体与外圈接触特性分析可知,滚动体在过故障期间与外圈接触间隙增大,接触力减小,并受到故障的阻碍作用,导致驱动力不足,因此其在过故障时与保持架兜孔两侧都会有碰撞,所以接触力有正有负[27]。对比图6中3种模型的滚动体过故障期间与保持架的接触力曲线可知,由于滚动体进出故障时有小幅度窜动,因此进出故障时与保持架有多个接触力峰值,且从图6中可以看出,故障塑性变形会导致碰撞次数增加,在图6中表现为,相对弹/弹模型、塑/塑模型和弹/塑模型滚动体进出故障时与保持架兜孔之间接触力波动次数更多,波动持续时间更长。

对比表6中不同模型进出故障时滚动体与保持架兜孔接触力平均值大小及变化范围可知,在滚动体进故障时,弹/弹模型滚动体和保持架接触力(绝对值)最大值80.97 N远小于塑/塑模型和弹/塑模型的该接触力(绝对值)最大值151.5 N及155.2 N,但是在出故障时则相反,因此弹/弹模型滚动体进、出故障时与保持架兜孔的接触力相差很大,其接触力平均值的(绝对值)差值为48 N;这是因为弹/弹模型滚动体进故障时故障边缘对滚动体作用力推动滚动体向前运动,出故障时弹/弹模型故障边缘对滚动体起阻碍作用;而塑/塑模型和弹/塑模型中故障处产生塑性变形,故障边缘形成一定的弧度,对滚动体切向作用力减小,因此这两种模型中滚动体进、出故障时与保持架兜孔接触力平均值的(绝对值)差值仅在1 N以内,由此可见故障塑性变形能够衰减滚动体与保持架兜孔前、后端接触力大小的不一致性;另外,根据上述故障处应变分析可知弹/塑模型故障处应变略大于塑/塑模型该处应变,这会导致滚动体与保持架接触力略微增大,因此表6中弹/塑模型滚动体在进、出故障时与保持架兜孔的接触力平均值的绝对值107.5 N、108.42 N要略大于塑/塑模型的96.85 N及96.35 N。

另外对比图6中3种模型的滚动体过故障时间可知,由于含塑性材料模型产生的塑性变形,导致滚动体与故障接触面积增大,滚动体提前进入故障,滞后出故障,在图6中表现为滚动体与保持架接触力最大值发生时间在进故障时提前,出故障时滞后,滚动体在故障中的时间延长。

表6 不同模型滚动体与保持架兜孔接触力统计值

4 结 论

本文以圆柱滚子轴承NU306为研究对象,在ANSYS LS-DYNA环境下建立了考虑滚动轴承故障塑性的3种动力学模型,分别为弹/弹模型、塑/塑模型、弹/塑模型,通过对3种轴承故障塑性动力学特性展开分析,得出主要结论如下。

对比实验与3种仿真模型分析结果,根据实验及理论的波形和外圈特征频率对比分析,验证了所建立轴承故障塑性动力学模型的有效性,进一步分析考虑滚动体产生塑性变形的塑/塑模型与不考虑滚动体产生塑性变形的弹/塑模型所得结果几乎一致,可知滚动体几乎不会产生塑性变形,因此在建模过程中可以仅考虑故障处的塑性变形。

轴承实际运行中故障产生后的局部接触引发应力集中,局部出现明显的塑性变形后,滚动体与故障处间隙增大,滚动体与故障处接触力显著减小,导致轴承中各滚动体受载荷不均,使得各滚动体运行状况不一致,轴承运行不稳定。

滚动体为塑性材料模型中,滚动体塑性变形后会出现塑性叠加现象,但是其变形量极小;故障塑性变形会引起滚动体与故障处接触面积增大,导致滚动体与保持架之间碰撞次数增多,故障边缘塑性变形成一定的弧度,使得滚动体与保持架兜孔前、后端接触力大小的不一致性减弱,接触力最大值发生的时间在进故障时会提前,出故障时会滞后。