多级齿轮传动裂纹故障频率分析与特征提取*

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

（新疆大学机械工程学院 乌鲁木齐，830047）

引言

齿轮传动系统具有高精度、高效率、高稳定性等优点而被广泛应用于航天航空、交通运输、能源装备等领域。由于齿轮传动系统结构日趋复杂，且工作环境相对恶劣，因此在其整个生命周期内易发生多种形式的故障。裂纹故障作为齿轮早期失效的主要形式，如何实现裂纹故障的判别与定位一直都是本研究领域中的重点。

针对此问题，国内外学者分别从故障机理及故障特征提取等方面开展了大量研究。Howard 等［1］利用有限元软件计算了裂纹故障下的齿轮副啮合刚度，并采用集中质量法建立了单级减速器的扭转振动模型。Wu 等［2］考虑了齿轮的横向振动建立了6自由度齿轮系统动力学模型，同时研究了含轮齿裂纹故障的传动系统的振动特性。在此基础上，万志国等［3］改进了啮合刚度的计算方法，并建立了含齿根裂纹传动系统的参数动力学模型。吴家腾等［4］采用了参数反求法对齿轮裂纹的时变啮合刚度进行求解，并将齿轮系统简化为包含2 个扭转刚度、4 个支撑刚度的动力学模型，求解了其振动响应，验证了反求法的有效性。Meng 等［5］用势能法计算了齿面剥落及裂纹故障下的齿轮啮合刚度，并分析了其故障特征。目前，时变啮合刚度求解方法多样且日益精确，但多数齿轮裂纹系统模型是采用集中质量法建立的，未有效考虑传动轴柔性，以致最后轴承位置振动响应的求解结果误差较大。为解决上述问题，常乐浩等［6］基于有限元法建立了单级斜齿轮传动系统动力学模型并求解其振动响应，结果表明，计入轴段柔性后齿轮副与轴承动态响应求解精度更高，更加贴合实际。

在故障特征提取方面，为解决传统信号分析方法缺乏自适应性而无法在强噪声、非线性的振动信号中提取到有效的故障特征信息，学者先后提出一系列信号处理办法并广泛应用于齿轮故障诊断中。于德介等［7］将经验模态分解（empirical mode decomposition，简称EMD）与自回归模型（auto-regressive，简称AR 模型）相结合，有效提取了齿轮故障的特征矢量。程军圣等［8］利用局部均值分解（local mean decomposition，简称LMD）与谱峭度（spectral kurtosis，简称SK）结合实现了对齿轮故障的识别。然而，EMD 和LMD 从本质上归属递归模式分解，存在模态混叠、端点效应、边界效应等问题［9-10］。近年来，Dragomiretskiy 等［11］提出了VMD，该方法不仅可避免EMD 与LMD 的不足，且对于噪声有良好的鲁棒性。Antoni［12］提出的快速谱峭度法可在强噪声中准确反映振动信号的瞬时冲击成分，以便快速选择最佳滤波器。

综上所述，为了研究含齿轮裂纹的多级齿轮传动系统故障特征，首先，采用有限元法求解两级齿轮副时变啮合刚度；其次，计入轴段柔性，建立二级齿轮传动系统有限元模型；最后，采集试验数据并用VMD-FSK 法对故障特征频率进行分析。

1 齿轮传动系统有限元模型

笔者以Spectra Quest 故障诊断综合试验台的二级圆柱直齿轮减速器传动系统为研究对象，传动系统由3 根转轴、2 对齿轮副以及3 对轴承组成。为提高轴承振动响应求解精度及故障振动信号的准确性，计入轴段柔性，建立传动系统的有限元模型如图1 所示，试验平台具体参数见表1，利用轴单元、齿轮啮合单元、轴承-轴单元将模型划分为32 个节点，29 个单元。轴承节点分别位于第5，13，14，22，23 和31 号节点，齿轮节点位于第7，16，20 和29 号节点，其中第7 号节点与第16 号节点组成第1 级齿轮啮合单元，第20 号节点与第29 号节点组成第2 级齿轮啮合单元，单元构成可参考文献［13］。

图1 二级齿轮传动系统有限元模型Fig.1 Finite element modeling of the two-stage gear transmission system

表1 仿真与试验齿轮参数Tab.1 Parameters of gear for simulation and test

2 裂纹故障下齿轮啮合单元动力学模型

在齿轮传动过程中，齿轮啮合刚度的周期性变化是造成齿轮箱振动、噪声的主要原因之一。齿根裂纹的产生会导致齿轮啮合刚度降低，从而导致齿轮箱异常振动。因此，建立含齿轮裂纹的动力学模型可为含轮齿裂纹的齿轮箱故障诊断提供理论基础。

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

建立如图2 所示的直齿圆柱齿轮6 自由度啮合单元动力学模型，考虑齿轮传动系统的时变啮合刚度Km（t）以及传递误差e（t）等对系统的影响。其中：vp，vg，wp，wg为齿轮副横向振动自由度；θp，θg为齿轮副扭转振动自由度。

图2 齿轮啮合单元动力学模型Fig.2 Dynamics model of meshing gear pair

将齿轮副各方向振动位移向啮合线方向投影得相对总变形量δ为

其中：rp，rg分别为主动轮与从动轮的基圆半径；α为压力角。

齿轮弹性啮合力可表示为

其中：km为齿轮副啮合刚度；cm为齿轮副啮合阻尼。

根据达朗贝尔原理，得到齿轮啮合单元的运动微分方程

其中：mp，mg分别为主、从动轮质量；Ip，Ig分别为主、从动轮转动惯量；fs为齿轮副的法向冲击力。

齿轮啮合单元的运动微分方程矩阵形式可表示为

其中：xm={vp，wp，θp，vg，wg，θg}T，为齿轮啮合单元位移向量；Fmo为齿轮啮合单元所受外部激励。

移项后，可得

其中：Fme为齿轮综合误差引起的内激励。

2.2 齿轮时变啮合刚度计算

建立直齿圆柱齿轮的三维接触模型，并采用有限元法对传动系统的齿轮啮合刚度进行求解。如图3 所示，在主动轮上施加负载扭矩，从动轮施加全约束，并采用罚函数法对接触进行定义。

图3 啮合刚度计算Fig.3 Calculation of meshing stiffness

在不考虑装配误差、加工误差以及齿侧间隙等因素的理想条件下，对齿轮副刚度进行求解。本模型采用六面体网格与四面体网格相结合的方式进行网格划分，采用六面体单元对轮齿与轮毂进行网格划分，其中非接触区域采用稀疏划分以适当减少计算量，在齿根位置处采用四面体单元进行过渡。

在齿轮啮合过程中常受到接触面之间的冲击影响，用罚函数法对接触进行定义，假设两接触点之间为线性弹簧连接［14］，若接触点之间相互渗透，根据胡克定律，两接触间相互作用力为

其中：k为接触刚度；μ为穿透深度。

接触刚度k为

其中：fo为罚函数系数，一般取0.1；M为接触单元体积模量；A为单元接触面积；V为接触单元体积。

对主动轮施加负载扭矩时，首先在局部圆柱坐标系中，对主动轮横向位移进行约束，再在齿轮内圈节点上施加节点力，以实现主动轮上负载扭矩的施加，扭矩计算公式为

其中：n为主动轮内圈节点数目；fm为周向节点力；rn为主动轮内圈半径。

在齿轮啮合线方向上产生的等效作用力为

其中：rb为主动轮基圆半径。

为避免局部变形的影响，提取主动轮内圈上所有节点位移，并取其平均值作为轮体内圈的扭转变形量δn，则主动轮的相对转角为

综合变形在啮合线方向上引起的等效变形为

齿轮啮合刚度计算公式可表示为

将主动轮某一轮齿的啮合周期划为10 等分，根据齿轮副啮合关系调整从动轮位置，获得齿轮副10种啮合状态，并重复以上刚度计算步骤得到10 种啮合状态对应的啮合刚度。此后，通过插值法拟合得到一个啮合周期的啮合刚度变化曲线如图3 所示。由图3 可知，整个啮合周期表现为双齿啮合-单齿啮合-双齿啮合。在单双齿交替时，啮合刚度出现阶跃式突变，造成传动系统呈现较为明显的刚度激励。

2.3 齿轮裂纹故障模型

轮齿在啮合过程中不断受到啮合力的冲击作用，在啮合面的齿根部位常常容易发生裂纹故障，构造裂纹如图4 所示。图4（a）中：d为齿宽；do为裂纹沿齿宽方向扩展长度；γ为裂纹角度；qo为裂纹深度。在裂纹尖端进一步加密网格以防止应力奇异，计算时重复啮合刚度计算步骤，最终拟合得到含齿根裂纹的齿轮啮合刚度曲线，如图4（b）所示。当第1 级齿轮副主动轮上发生裂纹故障时（γ=15ο），故障齿在啮合周期内刚度变化。与正常齿轮时变啮合刚度对比，可以看出，由于裂纹的存在导致了轮齿柔性增强，啮合刚度降低。随着裂纹深度的扩展，系统的啮合刚度也呈降低趋势，且裂纹对单齿啮合区的影响较双齿啮合区更为显著，单双齿交替突变幅值加大，从而在时域中表现为周期性突变。

图4 含齿根裂纹故障的轮齿模型Fig.4 Gear model with root crack fault

3 齿轮裂纹故障下传动系统模型总装

建立如文献［13］所示的级间相位关系，并综合考虑各单元的耦合关系与总刚度矩阵的储存带宽后，将各单元按节点编号进行集成，刚度总装示意图如图5 所示，整体矩阵对应系统节点自由度整体编号。去刚体位移后，总装矩阵包含95 个节点自由度，其中输入轴部分包含38 个节点自由度，中间轴部分包含27 个节点自由度，输出轴部分包含30 个节点自由度，空白处均为0 元素。

图5 含齿轮裂纹的二级齿轮传动系统刚度总装示意图Fig.5 Schematic diagram of integrated stiffness matrix

齿轮啮合单元作为连接轴系的耦合单元，当齿轮发生故障时，不仅导致齿轮所在节点的局部刚度降低，而且会使得耦合位置刚度值发生变化，从而实现故障在系统内的传递。

考虑含齿轮裂纹故障的二级齿轮传动系统的动力学微分方程表示为

其中：X为节点位移列阵；M，C，Kf均为95×95 矩阵，分别表示含齿轮裂纹故障系统的质量总装矩阵、系统的阻尼总装矩阵、系统的刚度总装矩阵；Po为系统外部激励；Fe为系统误差激励。

4 仿真分析

设输入转速w=500 r/min，负载扭矩T=100 N·m，可得输入轴转频ft1=8.3 Hz，中间轴转频ft2=3.3 Hz，输出轴转频ft3=0.96 Hz，第1 级齿轮副啮合频率fm1=299.8 Hz，第2 级齿轮副啮合频率fm2=95.7 Hz。通过求解其振动特征方程可得此时系统的固有频率，其中系统第1 阶固有频率为127.9 Hz。采用Newmark-β 法对传动系统动力学方程进行求解，提取输入轴右端的轴承振动信号如图6 所示。

图6 传动系统振动响应Fig.6 Vibration response of gear transmission system

图6（a）为无故障状态下传动系统的振动响应时域及频域图。由频域可知，系统振动响应以两级齿轮副啮合频率及其倍频为主，同时出现系统第1阶固有频率fn1。

图6（b）为在第1 级齿轮副主动轮存在裂纹时的振动响应时域及频域图，裂纹角度γ=15°，裂纹深度为2 mm。在时域中，出现明显的周期性冲击成分，冲击周期为1/ft1。分析其频域可知，在第1 级齿轮副啮合频率及其倍频附近产生大量边频，边频带间隔为ft1；而在第2 级齿轮副啮合频率及其倍频附近未发生明显变化，此时振动信号的载波频率为两级齿轮副啮合频率，调制频率为输入轴转频。

图6（c）为在第2 级齿轮副主动轮存在裂纹时的振动响应时域及频域图，裂纹参数同第1 级齿轮副主动轮。在时域中，冲击成分的冲击周期为1/ft2。分析其频域可知，边频成分集中在第2 级齿轮副啮合频率及其倍频附近，边频带间隔为ft2，而在第1 级齿轮副啮合频率及其倍频附近未发生明显变化，此时调制频率为中间轴转频。

图6（d）为两级齿轮副主动轮均存在裂纹时的振动响应时域及频域图。在时域中，冲击由两级齿轮副共同作用并相互传递，因此时域中出现的长冲击周期为1/ft1，短冲击周期为1/ft2。频谱图中，在输入轴转频及中间轴转频共同调制作用下，在第1 级齿轮副的倍频附近出现以ft1间隔的边频调制现象，同样在第2 级齿轮副的倍频附近也集中了大量的边频成分，其间隔为ft2。

基于此，可知转频fti（i=1，2，3）为传动系统齿轮裂纹的故障频率，当时域中出现以某转轴转动周期为冲击周期且在频域中出现某转轴的转频成分时，结合边频带中心可以实现故障齿轮的定位。

5 试验验证

本试验采用Spectra Quest故障诊断综合试验台，试验台如图7 所示，主传动链由二级定轴齿轮传动系统与单级行星齿轮传动系统构成。选用DEWESoft数据采集系统进行数据采集。在一级主动轮-轮齿齿根部位预制2 mm 裂纹，采样频率为20 480 Hz，转速为500 r/min，齿轮参数均与仿真一致。

图7 故障诊断综合试验平台Fig.7 Comprehensive test platform for fault diagnosis

图8 为试验台采集的原始振动信号，振动信号在时域上出现不规则冲击成分。在频域中可见振动信号边频带主要集中在齿轮啮合频率附近，但由于噪声等多因素影响，其故障特征无法明显体现。如图8（b）所示7 和20 倍的第1 级齿轮副啮合频率的边频带分布呈无规则分布，无法依此判断故障源，故需对振动信号进行处理。

图8 齿轮裂纹故障下振动响应Fig.8 Vibration response with gear crack fault

针对含故障的振动信号特点，首先选用VMD对原始信号进行分解，并采用相关系数法对分量进行筛选，最后利用FSK 与包络解调法对故障特征进行判别，具体处理流程如图9 所示。

图9 故障诊断流程图Fig.9 Flow chart of fault diagnosis

VMD 可以将振动信号f（t）分解为k个中心频率为ωk的模态分量，为此构建最优化变分模型如下

求解时，引入二次惩罚因子α与Lagrange 乘子λ（t）将原模型转化为非约束变分问题，得到扩展的Lagrange 表达式为

此后再利用乘子算法交替方向求解式（15），不断迭代，交替更新，λn+1如式（16）所示，以此搜索扩展Lagrange 的“鞍点”

此处对原始信号进行VMD 分解时，采用中心频率相近原则对分解层数K进行确定，此处K=9，惩罚因子α选择默认值2 000，τ=0，分解结果如图10 所示。

图10 VMD 处理结果Fig.10 Processing results by VMD

采用Pearson 相关系数对VMD 分量进行筛分，选取相关系数值最大的分量作为最优分量。其中第2分量相关度最高，因此选取第2分量继续执行下一步。

为实现谱峭度的快速计算，采用文献［15］中的FSK 对VMD 最优分量进行分析，得到快速峭度图如图11 所示。

图11 齿轮裂纹故障快速峭度图Fig.11 Fast kurtosis figure of gear crack fault

分析图11 可知，最大峭度值为24.7，对应在第4层分解信号。其滤波器参数如下：中心频率为7 360 Hz，带宽为640 Hz，依据该参数设置滤波器，完成滤波后结合包络分析得其信号包络谱，如图12所示。

图12 齿轮裂纹信号包络谱图Fig.12 The signal envelope spectrum for gear crack fault

由图12 可见包络谱清晰，故障特征频率及其倍频明显，故障频率为8.3 Hz，与输入轴转频一致，由此可以判断故障发生在第1 级齿轮副主动轮上。对比仿真分析，结论与仿真结果吻合。

6 结论

1）故障成分在系统振动响应的时域中表现出周期性冲击，冲击周期为故障齿轮所在轴转频的倒数。在频谱中，在频域中当某级齿轮副存在裂纹故障时，在该级齿轮副啮合频率及其倍频附近出现边频调制，其边频谱线间隔为故障齿所在轴的转频。

2）VMD-FSK 可有效实现多级齿轮副故障特征提取，试验结果与仿真得到的故障特征频率吻合。该方法能够较好地抵御噪声对信号的影响，可提高振动信号的信噪比。