计入传动轴柔性的齿轮故障系统建模及VMD-Lempel-Ziv损伤评估

申勇,章翔峰,姜宏,周建星,王成龙,乔帅,马铜伟

(新疆大学机械工程学院,830047,乌鲁木齐)

齿轮传动系统因具有高安全性、高稳定性等特点从而被广泛应用于电力、航天、能源装备等领域。然而,齿轮传动系统工作环境相对恶劣且受多源激励共同作用,故在其服役期间易发生不同程度故障,严重时常引发安全事故。因此,对其全生命周期劣化程度进行评估具有重要意义。

针对此课题,国内外学者分别从故障机理及故障诊断等方面开展了大量研究。Parker等利用半解析法求解了齿轮副二维有限元模型的啮合刚度及非线性特性[1-2]。Howard等借助有限元软件计算了裂纹故障下的齿轮副啮合刚度,并结合集中质量法建立了单级减速器的扭转振动模型[3]。在此基础上,Wu等考虑了齿轮的横向振动建立了6自由度齿轮系统动力学模型,研究了含轮齿裂纹故障的传动系统的振动特性[4]。张建宇等基于变截面阶梯悬臂梁假设建立含不同尺寸裂纹的轮齿模型,总结了轮齿劣化特性[5]。曲震等基于线弹性断裂力学研究了单对齿轮副的齿根裂纹轨迹及动态啮合过程的影响因素[6]。综上,在以往对于齿轮故障的研究,大多建立单对齿轮的集中质量法模型,未有效考虑传动轴柔性,以致得到的轴承位置振动响应及轴承间载荷分配误差较大。为此,Neriya较早将有限元法引入齿轮转子动力学领域[7]。在此基础上,常乐浩等基于有限元思想提出了一种齿轮箱通用建模方法,结果表明,计入轴段柔性后齿轮副与轴承动态响应求解精度更高,与现实情况更加贴合[8]。王胜男等分析了计入传动轴柔性后齿轮箱响应成分的变化,研究发现计入传动轴柔性后更能真实反映轴承间载荷分配关系[9]。

在以往的故障诊断研究中,大多构建其时频特征集以实现其故障诊断,并取得了一定成果。Dyer等统计了时域特征并构成了信号特征集,以此对轴承进行故障诊断[10]。而以往研究通常是将信号进行线性假设,易忽略其非线性特性。针对此类问题,Lempel等提出了Lempel-Ziv复杂度算法,该算法能较好地识别非线性系统信号的不规则成分,以此被广泛应用[11]。窦东阳等提出一种经验模态分解法(EMD)与Lempel-Ziv指标对轴承故障损伤程度进行判别[12]。张超等将局域均值分解法(LMD)与Lempel-Ziv指标对滚动轴承内外圈损伤程度进行判别[13]。EMD和LMD法从本质上归属递归模式分解,以致存在某些不可避免的问题,如模态混叠、端点效应、边界效应等问题[14-15]。近年来,Dragomiretskiy等提出了变分模态分解法(VMD),该方法不仅可避免EMD与LMD法的不足,且对于噪声有良好的鲁棒性,以被广泛应用于故障诊断领域[16]。

为了对齿轮损伤程度进行判别,本文采用有限元法求解了正常、裂纹故障、3mm断齿、9mm断齿4种状态下齿轮副时变啮合刚度;考虑传动轴柔性,建立二级齿轮传动系统有限元模型;提出一种VMD-Lempel-Ziv算法对试验信号的损伤程度进行判别。

1 二级齿轮传动系统模型简介

本文以SQI公司的故障诊断综合实验台中的二级直齿轮减速器传动系统为研究对象,如图1a所示。该系统由两对齿轮副、3根传动轴以及3组深沟球轴承组成。基于此,在建模过程中考虑传动轴柔性,采用有限元法建立该系统的动力学模型,模型单元划分如图1b所示。该模型共划分为32个节点、29个单元,主要由轴段单元及齿轮啮合单元构成,同时,考虑轴承约束作用建立轴-轴承耦合节点。其中,第1级齿轮啮合单元耦合作用于节点7与节点16,第2级齿轮啮合单元耦合作用于节点20与节点29,轴-轴承单元分别位于节点5、13、14、22、23、32。此外,齿轮故障位置均预设于第1级齿轮副的主动轮上。

(a)二级直齿轮减速器传动系统实物图

(b)二级齿轮传动系统有限元模型图1 二级齿轮传动系统模型

2 典型故障下传动系统动力学模型

2.1 齿轮副啮合单元动力学模型

在考虑齿轮时变刚度及综合误差等对啮合单元的影响下,构建6自由度啮合单元模型,具体齿轮副基本参数见表1。考虑齿轮啮合轮齿作用特点,将轮齿啮合视为黏弹性体作用,并用开尔文模型对其力学特性进行描述,即胡克体与牛顿体并联作用,齿轮啮合单元动力学模型如图2所示。

表1 齿轮副基本参数

图2 齿轮啮合单元动力学模型

图2中xp、yp、xg、yg表示齿轮副横向振动自由度,θp、θg为齿轮副扭转振动自由度,其中下标p与g分别代表主动轮与从动轮。齿轮副各自由度在啮合线方向投影得相对综合变形量为

δ=xpsinα+ypcosα-rpθp-xgsinα-

ygcosα-rgθg-e

(1)

式中:rp、rg分别为主、从动轮基圆半径;αo为压力角;e为综合传递误差。

结合达朗贝尔原理,得齿轮副振动微分运动方程为

(2)

式中:mp、mg分别为主、从动轮质量;cm为啮合阻尼;km(t)为齿轮啮合时变刚度;Ip、Ig分别为主、从动轮转动惯量;fs为齿轮副的法向冲击力。

由此,齿轮副啮合单元刚度矩阵Km可表示为

(3)

2.2 齿轮时变啮合刚度计算

齿轮啮合刚度周期性变化作为齿轮传动系统重要内部激励源,为探究其对系统振动影响,建立齿轮副三维接触有限元模型如图3所示。

图3 齿轮副三维接触有限元模型

在齿轮副啮合过程中轮齿接触之间受其法向冲击作用,将接触关系视为线弹性接触,并采用罚函数法对其进行定义[6]。

采用四面体网格与六面体网格相结合的方式对齿轮副进行网格划分。在轮齿及轮毂位置采用六面体网格进行划分,齿根位置选用四面体网格进行过渡。为减少计算量,在非接触区域进行稀疏划分。为保证计算可靠性对接触齿进行适当的网格加密。对从动轮内圈采用全约束,并在主动轮内圈节点上施加周向节点力,由此得到主动轮所受扭矩为

T=nfmrn

(4)

式中:n为主动轮内圈节点数;fm为周向节点力;rn为主动轮内圈半径。

在节点力作用下,主动轮存在扭转微位移,为避免局部变形影响,取主动轮内圈所有节点周向位移的平均值作为轮体节点变形量δn,因此主动轮相对转角可表示

φ=δn/rn

(5)

由此,可得在啮合线方向的等效作用力与综合变形量为

(6)

得到齿轮啮合刚度计算公式如下

(7)

图4 正常状态下齿轮副啮合刚度曲线

以第1级齿轮副为例,将主动轮某轮齿整个啮合周期进行等分,每两个状态转角相差1°,此后结合齿轮副传动比调整从动轮状态,以此得到该齿轮副10种啮合状态,并通过式(7)计算其啮合刚度。随后,采用插值法对散点数据进行拟合,以此方法得到两级齿轮副时变啮合刚度曲线如图4所示。由图可知,轮齿整个啮合周期表现为双齿啮合(啮入)-单齿啮合-双齿啮合(脱啮),从而整个传动过程啮合刚度呈阶跃式更替,以此成为齿轮传动过程中振动与噪声的主要来源之一。

2.3 齿轮典型故障模型

本文建立3种故障(裂纹、断齿宽度小于d/2、断齿宽度大于等于d/2)齿轮副模型,故障轮齿模型如图5所示。

(a)裂纹(b)断齿图5 故障轮齿模型

图5中d为齿宽;do为裂纹沿齿宽方向预置长度;γ为裂纹角度;qo为裂纹深度;qz为缺齿宽度。本文qo取2 mm,qz分别取3 mm、9 mm代表断齿宽度小于d/2或大于等于d/2。随后,针对以上轮齿缺陷程度,采用有限元法分别建立含不同故障的齿轮副三维接触模型,并结合2.2节方法对其啮合刚度进行计算,由此得到不同状态下的齿轮时变啮合刚度曲线如图6所示。

图6 不同状态下齿轮副时变啮合刚度曲线

由图6可知,由于故障的存在轮齿部分柔性增强,随着故障程度的加深,齿轮副啮合刚度呈不断下降的趋势。当断齿宽度大于等于d/2时,啮合刚度发生锐减,由于轮齿缺损较大,单对齿承载能力降低。此时,故障轮齿啮入状态下双齿啮合刚度趋向于正常状态下单齿啮合刚度数值。同时,正常齿与故障齿交替啮合,刚度发生异变,从而引起振动信号产生周期性冲击,冲击周期与故障齿交替速率相关。

2.4 轴承单元动力学模型

在齿轮传动系统中轴承主要起到支撑轴系及降低运动过程中的摩擦系数的作用。在轴承运转过程中,滚子依次通过受载区,其全过程为接触-接触形变最大-恢复变形,在该过程中轴承承载情况可分为奇压与偶压,轴承动力学模型如图7所示。

(a)奇压(b)偶压 图7 轴承动力学模型

奇压状态下,轴承刚度可表示为

(8)

式中:Fr为轴承的径向载荷;δr为该状态下轴承的径向位移,表达式如下

(9)

式中:Qmax为此时最大法向载荷;D为滚子直径;α为轴承的接触角。

将分配到滚子方位角±22.5°的载荷视为偶压状态最大载荷,并将该值代入式(8)(9)以得到轴承偶压状态下的刚度kbe。

将本模型中的3组深沟球轴承的刚度分布视为各向同性,则Kx(t)、Ky(t)分别代表水平方向与垂直方向轴承的时变刚度[17]

Kx,y(t)=Ko+Kasin(2πfbt+βo)

(10)

式中:Ko为轴承静刚度,Ko=(kbo+kbe)/2;Ka为轴承刚度波动幅值;fb为轴承通过频率;βo为轴承相位角。

3 模型总装

参照文献[17]考虑级间相位关系,结合各节点耦合关系及考虑刚度总装矩阵的储存带宽后,按各单元对应关系进行集成,含齿轮故障的二级齿轮传动系统总装示意图如图8所示。

图8 含齿轮故障的二级齿轮传动系统总装示意图

施加边界条件后,消除刚体位移的影响,得95×95规格的总体刚度矩阵。输入轴、中间轴及输出轴分别包括38、27、30个自由度。齿轮为轴间耦合的重要零部件,因此,当故障发生时,不仅影响齿轮节点位置刚度且易导致相啮合齿轮对应节点位置刚度发生变化。

总装后含齿轮故障系统的整体动力学微分方程表示为

(11)

式中:X表示节点位移列阵为{x1,y1,x2,y2,θ2,…,x31,y31,θ31,x32,y32,θ32};M、C、Kf分别表示含齿轮故障系统的质量总装矩阵、阻尼总装矩阵及刚度总装矩阵;Po为系统外部激励;Fe为系统误差激励。

4 系统振动响应及复杂度分析

4.1 时频分析

设置负载扭矩To为100 N·m,输入转速Ω为500 r/min,依据传动关系得输入轴转频ft1为8.3 Hz,中间轴转频ft2为3.3 Hz,输出轴转频ft3为0.96 Hz。结合齿轮参数,得啮合频率公式为

(12)

(13)

式中:fm1、fm2分别为第1级与第2级齿轮副啮合频率;Zp1、Zp2分别为第1级及第2级齿轮副主动轮齿数;Zg1为第1齿轮副从动轮齿数。由上得fm1=299.8 Hz,fm2=95.7 Hz。

为辨析系统响应频率成分以求解系统固有频率。求解时忽略阻尼对其振动响应的影响,因此,在外载荷F=0时,可得系统的自由振动方程

(14)

由此构造系统的振动特征方程

(15)

式中:wi为第i阶系统圆频率。

结合公式

(16)

将系统圆频率wi转化为系统固有频率,得到95阶系统固有频率,其中第1阶固有频率为fn1=127.9 Hz。

采用Newmark-β法对传动系统动力学方程式(11)进行求解,提取传动系统轴承振动响应如图9所示。

(a)正常

(b)裂纹

(c)断齿宽度3 mm

(d)断齿宽度9 mm图9 传动系统轴承振动响应

由图9可知,正常状态下系统响应成分主要包括两级齿轮副啮合频率及其倍频,同时在频率历程中亦出现系统第1阶固有频率。

在齿轮性能退化过程中,故障位置刚度发生异变,从而在系统响应时域图中出现周期性冲击,冲击周期为故障齿轮所在轴转频的倒数。在频域历程中,表现为故障齿轮啮合频率及其倍频为中心的频率调制,边频间隔为故障齿轮所在轴的转频。随着损伤程度的增大,时域历程中周期性冲击表现愈明显,响应幅值亦随之增高,在频域中,边频调制现象也越显著。

4.2 Lempel-Ziv复杂度指标分析

Lempel-Ziv指标是一种以符号序列表示时间序列的工具,以便对信号序列进行复杂度计算。本文结合提取的振动信号Z={z1,z2,z3,…,zN}进行复杂度计算,流程图如图10所示,具体步骤如下。

图10 Lempel-Ziv复杂度计算流程

(1)对振动信号进行二值粗粒化,即将信号进行[0,1]重构,得到符号序列S。为此,先对原信号序列取均值,并将信号中大于序列均值的点赋值为1,而小于均值的点赋值为0;

(2)初始化Po、Qo为空矩阵,令i=0,此时,复杂度C(i)=0;

(3)进入循环,令Pi={Pi-1si},Qi={Qi-1si},此后判断Pi-1是否包含Qi,若判断结果为“是”,复杂度C(i)保持不变,即C(i)=C(i-1)。而判断结果为“否”时,C(i)=C(i-1)+1,Qi={}。依此循环N次直至遍历序列S,通常N>3 600;

(4)复杂度CN归一化。对于二进制序列S,归一化复杂度计算如下

(17)

式中l=2。

不同状态下Lempel-Ziv复杂度变化仿真曲线如图11所示。其中,横坐标0代表正常状态,1表示裂纹状态,2为断齿3 mm,3代表断齿9 mm。

图11 不同状态下Lempel-Ziv复杂度变化仿真曲线

图11表明,系统出现不同程度故障时,其复杂度亦随之改变。故障状态下,系统出现转频调制,信号频率成分趋向复杂化。以致故障状态下的Lempel-Ziv指标整体高于正常状态。但随着损伤程度的递增,转频调制愈明显,频谱上啮合频率及其倍频附近边频带分布更加规律。因此,在齿轮退化过程中,轮齿发生早期故障时Lempel-Ziv指标增加最为明显,此后,随着轮齿损伤加剧,Lempel-Ziv指标逐渐减小。

5 试验验证

5.1 试验台简介及原始信号分析

本试验采用SQI公司的故障诊断综合试验台,如图12所示。该试验台主传动链由电机、两级平行轴齿轮变速箱以及单级行星齿轮变速箱构成,采用DEWESoft数据采集系统。并配备4种不同损伤状态(正常、裂纹、断齿宽度小于d/2、断齿宽度大于d/2)的第1级主动轮进行试验。采样频率设置为20 480 Hz、转速为500 r/min、负载为10 N·m,齿轮参数与仿真保持一致。

图12 故障诊断综合试验台

通过试验分析,得到原始试验信号,由于篇幅限制,本文以含断齿宽度大于d/2故障的振动信号为例进行分析,其原始试验信号的时频图如图13所示。

(a)时域响应

(b)频域响应图13 原始试验信号的时频图

由图13可知,其时域冲击成分显著,仿真分析与实验趋势吻合。在频域中,其幅值峰值主要集中在第1级齿轮啮合频率及其倍频附近,但边频带呈无规律分布,说明该信号受噪声等因素影响较大。此时,直接采用Lempel-Ziv指标对故障程度进行衡量,得到不同状态下未处理的试验信号Lempel-Ziv指标变化曲线如图14所示。

图14 不同状态下未处理的试验信号Lempel-Ziv 指标变化曲线

图14中,随着齿轮损伤程度的增加,但Lempel-Ziv值呈现增大—减小—增大的趋势,无法依此对齿轮损伤程度进行识别。

5.2 VMD-Lempel-Ziv算法及数据分析

针对5.1节信号特点,本文提出VMD-Lempel-Ziv算法对损伤进行评估,主要流程由以下几步构成:

(1)数据采集,设置VMD分解参数并对4种状态下的振动信号进行VMD分解;

(2)采用Person相关系数法,选取相关系数最大的IMF分量为最优分量;

(3)对最优IMF分量进行Lempel-Ziv计算,分析其损伤趋势。

VMD分解是以经典维纳滤波、Hilbert变换和频率混合为基础的变分问题求解过程[18]。VMD算法可将采集的振动信号f分解为K个以wk为中心频率的离散模态分量uk,从而构建变分模型如下

(18)

求解时,为将上述约束变分模型转化为非约束变分问题,引入二次惩罚因子α与拉格朗日乘子λ(t)构建扩展的拉格朗日表达式如下

L({uk},{wk},λ)=

(19)

VMD算法采用乘法算子交替方向法对式(19)进行求解,计算步骤如图15所示。

图15 VMD法的计算流程

(20)

(21)

(22)

本文依据文献[19]方法取分解层数K为5,惩罚因子α为2 000,τ为0。此处,以含断齿宽度大于d/2故障的振动信号为例,其分解结果如图16所示。

图16 VMD法的分解结果

此后,依据Person相关系数法选取最优IMF分量,并计算其Lempel-Ziv指标。依此,采用VMD-Lempel-Ziv法分别对4种典型故障振动信号进行计算,得到不同状态下试验信号Lempel-Ziv指标变化曲线如图17所示。

图17 不同状态下试验信号Lempel-Ziv指标变化曲线

分析图17可知,仿真计算结果与试验处理结果趋势基本一致,在故障早期系统振动信号Lempel-Ziv复杂度较正常情况呈增大趋势,而随着损伤程度的增加,振动信号的Lempel-Ziv值逐渐降低,与含齿轮故障的振动信号特征相符。同时,可发现经过VMD处理后的信号Lempel-Ziv复杂度整体小于图13中未经处理的Lempel-Ziv值,说明经VMD分解能消除部分噪声对振动特征的影响。

6 结 论

本文基于有限元法计算了齿轮副正常、裂纹故障、3 mm断齿、9 mm断齿4种状态下的时变啮合刚度,并计入传动轴柔性,建立了二级齿轮传动系统有限元动力学模型,最后采用VMD-Lempel-Ziv算法对试验信号进行判别,总结如下:

(1)齿轮故障主要致使振动信号发生转频调制,且随着故障程度的增加,其调制现象越明显,时域冲击幅值也随着增高。

(2)频率成分改变致使振动信号Lempel-Ziv复杂度亦随之变化,整个劣化过程的Lempel-Ziv复杂度呈先增后减的趋势,且在早期故障反应最为敏感。

(3)VMD-Lempel-Ziv算法能够有效的对齿轮副损伤程度进行判别,试验结果与仿真结果趋势基本吻合。同时,该方法可较好地降低噪声对信号的影响。