考虑齿面剥落的齿轮传动系统动力学建模

Dynamic modelling of gear transmission system considering tooth flank spalling

摘要

针对现有能量法无法准确描述齿轮存在齿面剥落的时变啮合刚度，造成齿轮传动系统故障的动力学建模方法不完善的问题，提出改进能量法。考虑不同类型的齿面剥落，建立一种考虑齿面剥落的时变啮合刚度公式。将所建立的齿面剥落时变啮合刚度公式引入六自由度传动系统，建立一个能够模拟齿面剥落的齿轮传动系统动力学模型。采用Newmark-β法进行求解，将数值模拟结果与试验数据进行对比分析，验证了所建立模型的准确性。基于该模型分析了不同齿面剥落关键参数对时变刚度的影响和对应的刚度变化规律及动力学响应的频率特性规律。研究结果表明，考虑齿面剥落的齿轮传动系统特征频率左右两侧会出现间距10Hz的次生频率峰值。

Abstract

Addressing the issue that the current energy method cannot accurately describe the time-varying meshing stiffness of the gears with tooth flank spalling，resulting in the imperfect dynamic modelling method for the fault of gear transmission system，an improved energy method was proposed.Considering various types of tooth flank spalling，a time-varying meshing stiffness formula that consi-dered the tooth flank spalling was established.The established time varying meshing stiffness was introduced into a six-degree-of-freedom transmission system，creating a dynamic model of the gear transmission system that could simulate the tooth flank spalling.The Newmark-β method was used for solving，the numerical simulation results were compared and analyzed with the experimental data to verify the accuracy of the established model.Based on this model，the influence of the main parameters of tooth flank spalling on the time-varying stiffness，the corresponding stiffness variation law and the frequency characteristics of the dynamic response were investigated.The results show that the secondary frequency peak with a spacing of 10 Hz occurs on both sides of the characteristic frequency of the gear transmission system with tooth flank spalling.

Graphical abstract

关键词

时变啮合刚度 / 齿面剥落 / 齿轮传动系统 / Newmark-β法 / 频率特性

Key words

time-varying meshing stiffness / tooth flank spalling / gear transmission system / Newmark-β method / frequency characteristics

齿轮传动系统在机械传动系统中应用十分广泛。根据齿轮传动特性，齿轮啮合时齿轮副的啮合刚度是动态变化的，齿轮的刚度变化是齿轮传动系统主要的内激励^［1］，因此时变啮合刚度是影响齿轮系统性能的重要因素之一^［2-3］。

研究齿轮时变啮合刚度的方法主要有能量法与有限元法两种。Dai等^［4］考虑负载变化，提出一种改进能量法计算行星齿轮中内齿轮的啮合刚度，并用试验进行了验证。Sun等^［5］基于薄片假设，建立了带齿形修形的直齿轮副时变啮合刚度修正模型，并在修正模型的基础上，研究了齿宽和扭矩对啮合刚度的影响。Huang等^［6］在切片法的基础上，提出了一种计算斜齿轮时变啮合刚度的新方法，采用“偏移和叠加”的思想避免重复计算每个切片齿轮的时变刚度，简化了计算流程。Wang等^［7］考虑轴向齿轮弯曲刚度、轴向齿轮扭转刚度和轴向齿轮基础刚度，提出一种改进的斜齿轮副时变啮合刚度计算方法。Zhao等^［8］考虑粗糙度对齿轮法向接触刚度的影响，建立了适合于齿轮副接触的分形接触模型，并基于此模型探究了粗糙度和加载静态传动误差对齿轮时变啮合刚度的影响。Marafona等^［9］总结了齿轮时变啮合刚度计算方法，将齿轮模型分为4种类型，分别是解析型、有限元型、混合型和近似解析型，为后续齿轮的研究指明了方向。

对于齿面剥落故障，学者们也逐渐通过改进模型精确度与提高计算精度等方法来得到更为精确的齿轮时变啮合刚度。Liang等^［10］研究了带点蚀外直齿轮的啮合刚度，模拟了不同形式点蚀对齿轮啮合刚度的影响。Ma等^［11］考虑扩展齿接触、双齿啮合区修正基体刚度、非线性接触刚度，分析齿轮剥落缺陷的影响，建立了齿轮剥落故障时变啮合刚度计算的解析模型。Luo等^［12］基于椭球几何结构，建立椭圆底面齿轮剥落模型并计算了啮合刚度。Lei等^［13］提出一种齿坑概率分布模型，用于计算带有齿面点蚀故障的啮合刚度。Luo等^［14］提出了一种基于缺陷率而不是特定的几何形状来模拟齿面剥落的方法。Wang等^［15］利用势能法，结合摩擦并考虑非线性Hertz接触刚度和修正圆角基础刚度，研究了带齿面剥落故障的齿轮啮合刚度。El等^［16］考虑现实的齿面剥落，研究了不同形状非均匀剥落和点蚀缺陷下的啮合刚度。Chen等^［17］研究了齿面磨损对齿轮时变啮合刚度的影响，并考虑了磨损的演变。上述文献的模型分析并不全面，有些没有深入考虑齿面剥落产生机理，有些则是没有找到齿面剥落关键参数对啮合刚度的影响，因此求解出的正常状态下及含有齿面剥落状态下的时变啮合刚度并不是十分精确。

现有对于齿轮传动系统的研究经历了由简单到复杂、由低自由度到高自由度的发展历程。He等^［18］分析了齿轮偏心对时变啮合刚度的影响，并采用改进的势能法建立新的齿轮啮合模型。在建立的十二自由度动力学模型的基础上，比较了有、无齿轮偏心时啮合刚度对齿轮系统的影响并定量研究了偏心对横向位移和扭转的响应及动态传动误差的影响。Xu等^［19］考虑了在困难工作条件下的旋转惯性、间隙和摩擦控制应用等影响因素，提出了一种分段非线性齿轮传动模型，并基于该模型对某重型机械臂的齿轮传动系统在启停和连续旋转条件下进行了实验研究。杨锐^［20］分别以直齿轮副和行星齿轮组为研究对象，开展了齿轮啮合动态激励的研究，在考虑了轮齿和轮体弹性变形的基础上，建立了健康齿轮和裂纹故障齿轮的啮合刚度模型，讨论了多自由度系统在刚度激励作用下的动态响应，并开展了针对齿轮系统故障机理和诊断关键技术的研究工作。Chen等^［21］考虑了单个齿轮的详细变形，包括轮齿柔性变形、齿轮体变形、局部轮齿接触挠度及由于相邻加载轮齿引起的轮齿变形，基于变形轮齿建立了正齿轮动力学模型，并与传统动力学模型进行了对比分析。但是，对于含有齿面剥落的齿轮传动系统动力学建模，现有研究存在齿面剥落机理不清晰、齿轮传动系统动力学响应分析不完善等不足。

本文针对时变啮合刚度求解精度不高与齿轮传动系统求解分析简单的问题，提出一种基于改进能量法的考虑齿面剥落故障的齿轮时变啮合刚度，并引入所建立的考虑弯扭耦合的六自由度齿轮传动系统动力学模型。对建立的模型用Newmark-β法进行求解，得到齿轮系统的动力学响应并进行时域频域分析。基于该模型与试验结果进行对比分析，验证了模型的准确性，得到了考虑齿面剥落的齿轮传动系统动力学响应特性。本文得到的结果可以为航空发动机齿轮箱故障诊断与寿命预测提供理论支撑。

1 齿面剥落时变啮合刚度

1.1 改进能量法时变啮合刚度

在实际生产与加工制造中，齿轮的齿根圆与基圆大部分情况下并不重合，而基础能量法没有考虑这一实际情况。根据机械原理，齿轮的齿根圆与基圆之间的关系为

R_{\mathrm{b}}=\frac{mz}{\mathrm{2}}{R}_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha =\frac{mz}{\mathrm{2}}-\left(h_{}^{*}+c^{\mathrm{*}}\right)m（1）

式中：R_{\mathrm{b}}为基圆半径；R_{\mathrm{f}}为齿根圆半径；m为模数；z为轮齿齿数；\alpha为齿轮压力角；h_{}^{*}为齿顶高系数；c^{\mathrm{*}}为顶隙系数。当取标准正齿轮齿顶高系数h_{}^{\mathrm{*}}=\mathrm{1}、顶隙系数c^{\mathrm{*}}=\mathrm{0.25}、压力角\alpha =\mathrm{20}°时，可求出当齿数约等于42时，齿根圆与基圆重合。

1.1.1 齿根圆半径小于基圆半径时的啮合刚度

当齿轮的齿根圆半径小于基圆半径时，即当轮齿齿数小于42时，轮齿简化计算模型如图1所示。

由图1可以看出，当齿根圆半径小于基圆半径时，之前的能量法计算公式只考虑超出基圆部分轮齿的刚度，并未考虑齿根圆与基圆中间部分的轮齿刚度。所以，需要结合齿轮齿形公式与能量法基本公式，对能量法计算公式进行修正，弯曲变形能修正为

U_{\mathrm{M}}={\int }_{\mathrm{0}}^d\frac{{\left[F_q\left(d-x\right)-F_{p}h\right]}^{\mathrm{2}}}{\mathrm{2}E{I}_x}\mathrm{d}x+{\int }_{\mathrm{0}}^{R_{\mathrm{b}}-R_f}\frac{{\left[F_q\left(d+x_{\mathrm{a}}\right)-F_{p}h\right]}^{\mathrm{2}}}{\mathrm{2}E{I}_{x_{\mathrm{a}}}}\mathrm{d}{x}_{\mathrm{a}}（2）

式中：x_{\mathrm{a}}为基圆往内的水平距离；I_{x_{\mathrm{a}}}为距基圆x_{\mathrm{a}}处的截面惯性矩。

因此改进能量法的弯曲刚度分量可表示为

\frac{\mathrm{1}}{k_{\mathrm{M}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}\frac{\mathrm{3}{\left\{\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}}\left[\mathrm{c}\mathrm{o}\mathrm{s}\alpha -\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]\right\}}^{\mathrm{2}}\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{\mathrm{2}EL{\left[\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}^{\mathrm{3}}}\mathrm{d}\alpha +{\int }_{\mathrm{0}}^{R_{\mathrm{b}}-R_{\mathrm{f}}}\frac{{\left[\left(d+x_{\mathrm{a}}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}}-h\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{1}}\right]}^{\mathrm{2}}}{E{I}_{x_{\mathrm{a}}}}\mathrm{d}{x}_{\mathrm{a}}（3）

类似地，剪切刚度、径向压缩刚度修正为

\frac{\mathrm{1}}{k_{\mathrm{S}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}\frac{\mathrm{1.2}\left(\mathrm{1}+v\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\alpha }_{\mathrm{1}}\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{EL\left[\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}\mathrm{d}\alpha +{\int }_{\mathrm{0}}^{R_{\mathrm{b}}-R_{\mathrm{f}}}\frac{{\left(\mathrm{1.2}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}}\right)}^{\mathrm{2}}}{G{A}_{x_{\mathrm{a}}}}\mathrm{d}{x}_{\mathrm{a}}（4）

\frac{\mathrm{1}}{k_{\mathrm{N}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{2}}}\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\alpha }_{\mathrm{1}}\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{EL\left[\left({\alpha }_{\mathrm{2}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}\mathrm{d}\alpha +{\int }_{\mathrm{0}}^{R_{\mathrm{b}}-R_{\mathrm{f}}}\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\alpha }_{\mathrm{1}}}{E{A}_{x_{\mathrm{a}}}}\mathrm{d}{x}_{\mathrm{a}}（5）

1.1.2 齿根圆半径大于基圆半径时的啮合刚度

当齿根圆半径大于基圆半径，即轮齿齿数大于42时，轮齿简化计算模型如图2所示。

由图2可以看出，基础能量法考虑了从基圆延伸出的轮齿部分，但是当齿数大于42时，真正的轮齿部分是从齿根圆延伸出的，所以基础能量法公式同样需要进一步修正。

结合图2，各刚度分量表达式可修正为

\frac{\mathrm{1}}{k_{\mathrm{M}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{3}}}\frac{\mathrm{3}{\left\{\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{1}}\left[\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]\right\}}^{\mathrm{2}}\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{\mathrm{2}EL{\left[\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}^{\mathrm{3}}}\mathrm{d}\alpha（6）

\frac{\mathrm{1}}{k_{\mathrm{S}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{3}}}\frac{\mathrm{1.2}\left(\mathrm{1}+v\right)\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\alpha }_{\mathrm{1}}\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{EL\left[\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}\mathrm{d}\alpha（7）

\frac{\mathrm{1}}{k_{\mathrm{N}}}={\int }_{-{\alpha }_{\mathrm{1}}}^{{\alpha }_{\mathrm{3}}}\frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\alpha }_{\mathrm{1}}\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha }{\mathrm{2}EL\left[\left({\alpha }_{\mathrm{3}}-\alpha \right)\mathrm{c}\mathrm{o}\mathrm{s}\alpha +\mathrm{s}\mathrm{i}\mathrm{n}\alpha \right]}\mathrm{d}\alpha（8）

式中：{\alpha }_{\mathrm{3}}为齿轮齿根圆上轮齿所占圆心角的一半。

利用改进能量法所建立的齿轮时变啮合刚度计算方法，可以计算出齿轮啮合时的时变啮合刚度变化曲线。图3是应用改进能量法计算结果绘制的时变啮合刚度图像。

图4是包含2个双齿啮合区与2个单齿啮合区的时变啮合刚度放大图像。由图3和图4可以看出与基础能量法一样的刚度变化规律，但改进能量法的单齿啮合区域比基础能量法计算得到的时变啮合刚度曲线更凸出一些。

1.2 改进能量法与基础能量法的对比分析

为了验证改进能量法计算得出的齿轮时变啮合刚度的准确性与合理性，将基础能量法与改进能量法得出的时变啮合刚度进行比较。另外，由于石川公式法计算得出的齿轮时变啮合刚度已被确定为国际标准，所以也可作为对照基准纳入比较（石川公式法的详细介绍参见文献［22］）。

齿数等于42时，基础能量法与改进能量法得出的结果相同，所以只对比分析齿数小于42和大于42时的情况，分别取齿数为20和60。由图5可知，基础能量法在单齿啮合区比改进能量法计算结果更接近石川公式法，而改进能量法的计算结果则在双齿啮合区更接近石川公式法。对此进行定量分析，取单齿啮合区和双齿啮合区的刚度平均值并分析相对误差作为评判标准。

表1为不同刚度求解方法对比表。由表1可以看出，基础能量法与改进能量法在单齿啮合区和双齿啮合区的误差各不相同。齿数为20时，单齿啮合区基础能量法比改进能量法更接近石川公式法，双齿啮合区改进能量法比基础能量法更接近石川公式法；齿数为60时，单齿啮合区与双齿啮合区改进能量法都比基础能量法效果更好。取平均相对误差后可以看出，改进能量法比基础能量法更接近石川公式法，验证了改进能量法求解齿轮啮合刚度的有效性与准确性。

1.3 基于改进能量法的齿面剥落时变啮合刚度

1.3.1 齿轮齿面剥落时变啮合刚度机理

由于齿轮啮合时会产生一个很强的瞬时应力，所以在轮齿接触位置往往会造成齿面剥落。齿面剥落常见于轮齿中部，主要是由于轮齿长时间受到循环剪应力载荷，导致齿轮承受应力超过轮齿材料疲劳极限，齿轮产生微小裂纹。在冲击力和滑油等的作用下，微裂纹逐步扩大最终造成小块金属从齿面上脱落，形成齿面剥落或点蚀。

图6为齿面剥落区块示意图。如图6所示，由于齿面剥落区块尺寸较小，轮齿整体的结构未发生较大变化，所以对轮齿的弯曲刚度、径向压缩刚度和剪切刚度的影响均较小。但是由于齿面产生了剥落区块，齿轮齿廓发生了变化，表面接触特性也发生了改变，进而影响了赫兹接触刚度。因此，齿轮总刚度仍为上述各刚度分量的总和，其中赫兹接触刚度变为考虑剥落故障的赫兹接触刚度，其余各刚度分量保持不变。

假设矩形剥落区块的深度为0.1 mm，长度为L_sp，宽度为H_sp，则啮合齿轮的有效线长度变为

L_{\mathrm{s}}=L-L_{\mathrm{s}\mathrm{p}}（9）

相应的，赫兹接触刚度变为

\frac{\mathrm{1}}{k_{\mathrm{h}\_\mathrm{s}\mathrm{p}}}=\frac{\mathrm{4}\left(\mathrm{1}-{\nu }^{\mathrm{2}}\right)}{\mathrm{\pi }E{L}_{\mathrm{s}\mathrm{p}}}（10）

1.3.2 齿面剥落对时变啮合刚度的影响分析

由上述推导及图6可以分析得出，决定齿面剥落的关键参数在于剥落区块的长度L_sp及宽度H_sp，下面分别分析其对啮合刚度的影响。

（1）剥落区块长度对齿轮时变啮合刚度的影响

固定剥落区块宽度H_sp为5 mm，假设其落在啮合齿轮的某双啮合区间中，分别取5种剥落区块长度L_sp，绘制出相应的时变啮合刚度图像，如图7所示。

分析图7可知，齿面剥落长度对啮合刚度影响是在局部改变啮合刚度的大小。在齿面剥落宽度相等的情况下，剥落区块越长，刚度减小得越多，并且从图7可以看出，随着齿轮齿数增多，齿轮刚度减小的幅度逐渐减小。由于剥落位置的宽度不变，所以齿轮齿数越大，其对单个轮齿刚度的影响就越小。

（2）剥落区块宽度对齿轮时变啮合刚度的影响

固定剥落区块长度L_{\mathrm{s}\mathrm{p}}为12 mm，分别取5种剥落区块宽度H_{\mathrm{s}\mathrm{p}}，绘制出相应的齿轮时变啮合刚度图像，如图8所示。

由图8可以看出，齿面剥落宽度影响齿面啮合刚度下降的区间宽度，并且由于齿面剥落区块较小，对整体啮合刚度的影响不明显，仅在剥落区块所在区间上导致时变啮合刚度减小。随着齿面剥落宽度的增大，齿轮啮合刚度下降区间的宽度也增大。齿面剥落对于齿轮时变啮合刚度的影响相较其他几种故障类型更小，这是因为齿面剥落产生的点蚀区块尺寸较小，对轮齿结构的影响也较小，但是会造成齿面的应力集中，时间长了也会造成重大危害。

2 考虑齿面剥落的齿轮传动系统动力学建模

2.1 齿轮传动系统动力学建模

如图9所示，将一对啮合齿轮副离散为两个集中质量圆盘，不考虑支承轴承的影响，系统包含主动轮与从动轮的x、y方向平动自由度与主动轮和从动轮的两个转动自由度，主动轮与从动轮之间由啮合刚度与啮合阻尼连接，最终建立一个齿轮传动系统六自由度扭转振动模型。图9中下角标1代表主动轮，下角标2代表从动轮，k代表刚度，c代表阻尼。

对主动轮和从动轮进行运动学与动力学分析，列动力学微分方程，得到

\left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{x}}_{\mathrm{1}}-c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}+c_{\mathrm{1}}{\dot{x}}_{\mathrm{1}}-k_{\mathrm{12}}{\delta }_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}+k_{\mathrm{1}}{x}_{\mathrm{1}}=F_{x\mathrm{1}}\\ m_{\mathrm{1}}{\ddot{y}}_{\mathrm{1}}-c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}+c_{\mathrm{1}}{\dot{y}}_{\mathrm{1}}-k_{\mathrm{12}}{\delta }_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}+k_{\mathrm{1}}{y}_{\mathrm{1}}=F_{y\mathrm{1}}\\ I_{\mathrm{1}}{\ddot{\theta }}_{\mathrm{1}}+c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}{r}_{b\mathrm{1}}+k_{\mathrm{12}}{\delta }_{\mathrm{12}}{r}_{b\mathrm{1}}=T_{\mathrm{i}\mathrm{n}}\\ m_{\mathrm{2}}{\ddot{x}}_{\mathrm{2}}+c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}+c_{\mathrm{2}}{\dot{x}}_{\mathrm{2}}+k_{\mathrm{12}}{\delta }_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}+k_{\mathrm{2}}{x}_{\mathrm{2}}=F_{x\mathrm{2}}\\ m_{\mathrm{2}}{\ddot{y}}_{\mathrm{2}}+c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}+c_{\mathrm{2}}{\dot{y}}_{\mathrm{2}}+k_{\mathrm{12}}{\delta }_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}+k_{\mathrm{2}}{y}_{\mathrm{2}}=F_{y\mathrm{2}}\\ I_{\mathrm{2}}{\ddot{\theta }}_{\mathrm{2}}+c_{\mathrm{12}}{\dot{\delta }}_{\mathrm{12}}{r}_{b\mathrm{2}}+k_{\mathrm{12}}{\delta }_{\mathrm{12}}{r}_{b\mathrm{2}}=T_{\mathrm{o}\mathrm{u}\mathrm{t}}\end{array}\right.（11）

式中：x、y分别为两个齿轮x、y方向上的位移；\theta为转角；k_{\mathrm{12}}为齿轮时变啮合刚度；c_{\mathrm{12}}为齿轮啮合阻尼；F_{x\mathrm{1}}与F_{y\mathrm{1}}分别为主动轮x方向与y方向的外激励；F_{x\mathrm{2}}与F_{y\mathrm{2}}分别为从动轮对应方向的外激励；T_{\mathrm{i}\mathrm{n}}与T_{\mathrm{o}\mathrm{u}\mathrm{t}}分别为输入扭矩与输出扭矩，其余参数详见下文分析表述。

此外，上述动力学方程中

{\delta }_{\mathrm{12}}=\left(x_{\mathrm{1}}-x_{\mathrm{2}}\right)\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}-\left(y_{\mathrm{1}}-y_{\mathrm{2}}\right)\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}+{\theta }_{\mathrm{1}}{r}_{b\mathrm{1}}+{\theta }_{\mathrm{2}}{r}_{b\mathrm{2}}+e\left(t\right)（12）

式中：{\alpha }_{\mathrm{0}}为齿轮啮合压力角；r_{b\mathrm{1}}与r_{b\mathrm{2}}分别为主动轮和从动轮的基圆半径；e(t)为齿轮的静态传递误差。下面分别分析动力学方程中的质量矩阵、阻尼矩阵和刚度矩阵。

（1）质量矩阵

系统的质量矩阵可表示为主动轮和从动轮的两个涉及平动的质量与涉及转动的惯量，矩阵表达式为

\left\{\begin{array}{l}\mathit{M}=\left[\begin{array}{cc}{M}_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& M_{\mathrm{2}}\end{array}\right]\\ {\mathit{M}}_{\mathrm{1}}=\left[\begin{array}{ccc}{m}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& m_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& I_{\mathrm{1}}\end{array}\right]\\ {\mathit{M}}_{\mathrm{2}}=\left[\begin{array}{ccc}{m}_{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& m_{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& I_{\mathrm{2}}\end{array}\right]\end{array}\right.（13）

式中：m为质量；I为转动惯量；下角标1、2分别为主动轮与从动轮。

（2）阻尼矩阵

本文采用文献［23］中给出的齿轮啮合黏性阻尼的经验公式

c_{\mathrm{12}}=\mathrm{2}{\xi }_m\sqrt[]{k_m\frac{r_{b\mathrm{1}}^{\mathrm{2}}{r}_{b\mathrm{2}}^{\mathrm{2}}{I}_{\mathrm{1}}{I}_{\mathrm{2}}}{I_{\mathrm{1}}{r}_{b\mathrm{1}}^{\mathrm{2}}+I_{\mathrm{2}}{r}_{b\mathrm{2}}^{\mathrm{2}}}}（14）

式中：k_m为啮合刚度；{\xi }_m为啮合阻尼比，文献［23］给出的参考值为0.03~0.17，本文取啮合阻尼比为0.1。带入相关参数，啮合刚度取双齿啮合与单齿啮合的平均刚度，计算可得相关阻尼参数。

对式（11）进行整理，可得到此动力学系统的阻尼矩阵\mathit{C}，即

\mathit{C}=\left[\begin{array}{llllll}{c}_{\mathrm{12}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\alpha +c_{\mathrm{1}}& c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & -c_{\mathrm{12}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\alpha & -c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \\ c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\alpha +c_{\mathrm{1}}& -c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & -c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -c_{\mathrm{12}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\alpha & -c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \\ -c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & -c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}^{\mathrm{2}}& c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}{r}_{b\mathrm{2}}\\ -c_{\mathrm{12}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\alpha & -c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & c_{\mathrm{12}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\alpha +c_{\mathrm{2}}& c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \\ -c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & -c_{\mathrm{12}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}+c_{\mathrm{2}}& c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha \\ -c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & -c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{1}}{r}_{b\mathrm{2}}& c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\alpha & c_{\mathrm{12}}{r}_{b\mathrm{2}}^{\mathrm{2}}\end{array}\right]（15）

由式（15）可知，该系统阻尼矩阵为对称阵，包含的关键参数为2个齿轮的支承阻尼、齿轮副的啮合阻尼及2个齿轮的基圆半径和齿轮的压力角。

（3）刚度矩阵

在此动力学模型中，刚度由齿轮的支承刚度与齿轮副啮合刚度组成。其中啮合刚度在第1节中已经给出，齿轮支承刚度采用参考文献［24］中给出的值。刚度矩阵的表述与阻尼矩阵完全相同，仅其中的阻尼参数被替换为相对应的刚度参数。本文用到的所有齿轮质量阻尼及刚度参数如表2所示。

2.2 齿面剥落传动系统动力学建模与响应特性

本文采用Newmark-β法对系统进行求解。选取Newmark-β法的计算参数\varsigma为0.5，β为0.25，即采用平均加速度法来求解动力学模型。

齿轮传动系统动力学方程可表示为

\mathit{M}\ddot{q}+\mathit{C}\dot{q}+\mathit{K}\mathit{q}=\mathit{F}（16）

式中：\mathit{M}为系统的质量矩阵；\mathit{C}为系统的阻尼矩阵；\mathit{K}为系统的刚度矩阵；\mathit{F}为系统的外部激励矩阵；\mathit{q}是系统的广义坐标矩阵。对此动力学方程应用Newmark-β法求解，可得系统的广义坐标响应矩阵。齿轮啮合频率f_m为

f_m=z_{\mathrm{1}}\cdot n_{\mathrm{1}}/\mathrm{60}=z_{\mathrm{2}}\cdot n_{\mathrm{2}}/\mathrm{60}（17）

式中：z_{\mathrm{1}}和z_{\mathrm{2}}为主动轮和从动轮的齿数；n_{\mathrm{1}}为输入电机转速；n_{\mathrm{2}}为输出负载轴转速。

2.2.1 正常齿轮传动系统动力学响应分析

设置输入电机转速为1 200 r/min，当主动轮与从动轮齿数为20时，齿轮啮合频率为400 Hz；当主动轮与从动轮齿数为42时，齿轮啮合频率为840 Hz。由本文中模型计算得到的时域频域图如图10所示。

设置输入电机转速为1 800 r/min，当主动轮与从动轮齿数为20时，齿轮啮合频率为600 Hz；当主动轮与从动轮齿数为42时，齿轮啮合频率为1 260 Hz。由本文中模型计算得到的时域频域图如图11所示。

2.2.2 齿面剥落齿轮传动系统动力学响应分析

将正常情况下的齿轮时变啮合刚度替换为齿面剥落的齿轮时变啮合刚度，即可得到带有齿面剥落故障的齿轮传动系统动力学模型。

1.3节得到关于齿面剥落有剥落区块长度L_{\mathrm{s}\mathrm{p}}与剥落区块宽度H_{\mathrm{s}\mathrm{p}}，分别对应数值为8 mm与5 mm，取齿轮齿数为32，分别取齿轮转速为1 200 r/min和1 800 r/min，即转频为20 Hz和30 Hz，得到如图12的齿面为32时不同转速下的时频域图。

由图12可以看出，传动系统的时域呈周期性变化，且在周期开始阶段振动加速度最大；在频域图中，除了出现齿轮啮合特征频率，还在其左右两侧出现了次生频率，其与峰值频率的间距为10 Hz。

2.3 试验验证

本文采用东南大学齿轮数据集进行对比验证分析，关于齿轮箱的详细参数见文献［25］。齿轮箱和振动试验台如图13所示。

该试验台由电机、电机控制器、行星齿轮箱、减速齿轮箱、刹车和刹车控制器组成。试验台可进行轴承与齿轮的故障数据采集，也可以进行整体齿轮-轴承系统动力学数据采集。试验的采样频率为5 120 Hz，一共有8个通道信号，具体如下：第1列为输入电机振动信号；第2、3、4列分别对应行星齿轮x、y和z 3个方向的振动信号；第5列对应着电机负载信号；第6、7、8列分别对应着减速器x、y和z 3个方向的振动信号。

下面分别对含有齿面剥落的齿轮转频为20 Hz和30 Hz 2种工况，选取七通道传感器数据，即减速器y方向的振动信号进行分析。减速器位置对应的齿轮齿数为32，得到对应的时频域图如图14所示。

由图14可以看出，当啮合齿轮发生齿面剥落故障时，齿轮振动信号时域图出现周期性变化规律，并且周期开始时会出现强振动信号。齿轮频域图则出现齿轮啮合特征频率，并且在其左右两侧均出现了次生频率峰值，间距为10 Hz。由东南大学试验台得出的含齿面剥落的齿轮传动系统动力学响应与理论模型计算的齿轮系统动力学响应规律一致，验证了本文提出的含有齿面剥落的齿轮传动系统动力学模型的准确性与合理性。

3 结论

本文针对现有的含有齿面剥落的齿轮啮合时变啮合刚度存在的求解方法不充足及含有齿面剥落的齿轮传动系统动力学建模研究不足等问题，提出一套齿面剥落故障的齿轮时变啮合刚度求解公式，分析了不同参数齿面剥落对啮合刚度的影响，并将其引入六自由度扭转振动动力学模型。采用Newmark-β法求解出其动力学响应，并对其进行傅里叶变换并进行频谱分析。本文最终得到如下的结论：

（1）建立了一个六自由度扭转振动齿轮传动系统动力学模型，引入含齿面剥落的齿轮时变啮合刚度并分析其对模型动力学响应的影响，与东南大学试验台试验结果对比分析，验证了本文模型的合理性与准确性。

（2）分析了齿面剥落的不同关键参数对齿轮时变啮合刚度的影响规律。齿面剥落长度对啮合刚度影响的方式是局部改变啮合刚度的大小，在齿面剥落宽度相等的情况下剥落区块越长，刚度减少得越多；齿面剥落宽度影响的是齿面啮合刚度下降的区间宽度，并且由于齿面剥落区块较小，对整体啮合刚度的影响不明显，仅在剥落区块所在区间上导致时变啮合刚度的减少。

（3）当齿轮正常时，齿轮振动信号时域图出现较为明显的周期性，而齿轮振动信号经过频谱分析后会出现齿轮啮合特征频率，其大小为齿轮齿数乘以齿轮轴转频。而当齿轮发生齿面剥落时，齿轮振动信号时域图会出现较为明显的周期性，并且振动信号加速度最大的部分出现在啮合周期前部。经过傅里叶变换得出频域图以后，除了齿轮啮合频率以外，啮合特征频率左右侧均出现了次生频率，并且相距为10 Hz。

参考文献

[1] Xue S， Howard I， Wang C S，et al.Dynamic mode-lling of the gear system under non-stationary conditions using the iterative convergence of the tooth mesh stiffness［J］.Engineering Failure Analysis，2022，131：105908.

[2] Sun Z， Chen S Y， Tao X，et al.Research on dynamic characteristics of gear system considering the influence of temperature on material performance［J］.Journal of Vibration and Control，2022，28（13/14）：1792-1803.

[3] Wang J， Zhang J.Effects of random interval para-meters on spur gear vibration［J］.Journal of Vibration and Control，2021，27（19/20）：2332-2344.

[4] Dai H， Long X H， Chen F，et al.An improved analytical model for gear mesh stiffness calculation［J］.Mechanism and Machine Theory，2021，159（1）：104262.

[5] Sun Y N， Ma H， Huang Y F，et al.A revised time-varying mesh stiffness model of spur gear pairs with tooth modifications［J］.Mechanism and Machine Theory，2018，129（8）：261-278.

[6] Huang Y F， Chen K K， Ma H，et al.Deformation and meshing stiffness analysis of cracked helical gear pairs［J］.Engineering Failure Analysis，2019，95：30-46.

[7] Wang Q B， Zhao B， Fu Y，et al.An improved time-varying mesh stiffness model for helical gear pairs considering axial mesh force component［J］.Mechanical Systems and Signal Processing，2018，106（1）：413-429.

[8] Zhao Z F， Han H Z， Wang P F，et al.An improved model for meshing characteristics analysis of spur gears considering fractal surface contact and friction［J］.Mechanism and Machine Theory，2021，158（1）：104219.

[9] Marafona J D M， Marques P M T， Martins R C，et al.Mesh stiffness models for cylindrical gears：a detailed review［J］.Mechanism and Machine Theory，2021，166：104472.

[10] Liang X H， Zhang H S， Liu L B，et al.The influence of tooth pitting on the mesh stiffness of a pair of external spur gears［J］.Mechanism and Machine Theory，2016，106（8）：1-15.

[11] Ma H， Li Z W， Feng M J，et al.Time-varying mesh stiffness calculation of spur gears with spalling defect［J］.Engineering Failure Analysis，2016，66（4）：166-176.

[12] Luo Y， Baddour N， Han G S，et al.Evaluation of the time-varying mesh stiffness for gears with tooth spalls with curved-bottom features［J］.Engineering Failure Analysis，2018，92（6）：430-442.

[13] Lei Y G， Liu Z Y， Wang D L，et al.A probability distribution model of tooth pits for evaluating time-varying mesh stiffness of pitting gears［J］.Mechanical Systems and Signal Processing，2018，106（1）：355-366.

[14] Luo Y， Natalie B， Liang M.A shape-independent approach to modelling gear tooth spalls for time varying mesh stiffness a evaluation of a spur gear pair［J］.Mechanical Systems and Signal Proce-ssing，2009，120：836-852.

[15] Wang S Y， Zhu R P.An improved mesh stiffness calculation model of spur gear pair under mixed EHL friction with spalling effect［J］.Vibroengineering Procedia，2020，33（2）：176-181.

[16] Elyousfi B， Soualhi A， Medjaher K，et al.New app-roach for gear mesh stiffness evaluation of spur gears with surface defects［J］.Engineering Failure Analysis，2020，116（4）：104740.

[17] Chen W， Lei Y L， Fu Y，et al.A study of effects of tooth surface wear on time-varying mesh stiffness of external spur gear considering wear evolution process［J］.Mechanism and Machine Theory，2021，155（1）：104055.

[18] He X Z， Zhou X Q， Xue Z，et al.Effects of gear eccentricity on time-varying mesh stiffness and dynamic behavior of a two-stage gear system［J］.Journal of Mechanical Science and Technology，2019，33（3）：1019-1032.

[19] Xu Z， Deng H， Zhang Y.Piecewise nonlinear dynamic modeling for gear transmissions with rotary inertia and backlash［J］.IEEE Access，2019，7：176495-176503.

[20] 杨锐.齿轮非线性动态特性分析及裂纹故障诊断研究［D］.南宁：广西大学，2021.

[21] Chen Z G， Ning J Y， Wang K Y，et al.An improved dynamic model of spur gear transmission considering coupling effect between gear neighbo-ring teeth［J］.Nonlinear Dynamics，2021，106（1）：339-357.

[22] Wang J G， He G Y， Zhang J，et al.Nonlinear dynamics analysis of the spur gear system for railway locomotive［J］.Mechanical Systems and Signal Processing，2017，85（8）：41-55.

[23] Lin H H， Huston R L， Coy J J.On dynamic loads in parallel shaft transmissions：part I-modelling and analysis［J］.Journal of Mechanisms，Transmi-ssions，and Automation in Design，1988，110（2）：221-225.

[24] Fang B， Zhang J H， Yan K，et al.A comprehensive study on the speed-varying stiffness of ball bea-ring under different load conditions［J］.Mechanism and Machine Theory，2019，136：1-13.

[25] 陈超.基于迁移学习的旋转机械故障诊断方法研究［D］.南京：东南大学，2020.

基金资助

国家自然科学基金(12172231)

辽宁省兴辽英才计划项目(XLYC2203042)

沈阳市中青年科技创新人才支持计划项目(RC220439)