基于齿轮振动和非线性接触的动态响应理论建模

付 饶,李静宇,熊奉奎

(1.绵阳职业技术学院 机电工程系，四川 绵阳 621000； 2.沈阳农业大学,辽宁 沈阳 110000；3.陕西科技大学 机电工程学院，陕西 西安 710016 )

0 引 言

近年来，齿轮传动系统动态响应和振动分析，已经成为一个热门研究方向，人们对此进行了较为广泛研究。齿轮传动系统是一种自激系统，其中固有因素和外部因素共同激发了传动系统的振动，再加上具有自激机制的周期性接触损失，便会导致包含周期性振动和混沌响应的高度非线性响应的出现[1]。通常认为，通过改善齿轮的制造误差、接触齿对的数量及接触位置，以改变啮合刚度在啮合循环中的变化是优化传动系统响应和振动特性的最佳手段[2]。

基于FEA引入的计算模型，Parker等[3]，首先利用有限元分析来获得预先计算的齿轮啮合刚度，再通过引入一种半分析方法进行分析，使计算效率得到了显着提高。此种半分析方法使用近场解析方案替代了FEA的远场方案，使模型准确再现了Kahraman和Blankenship[4]的实验结果，但计算精度相对降低。随后Byrtus和Zeman[5]在此基础上提出了包含轴承和轴承座的齿轮箱分析模型，并通过结构优化和降低自由度，得到了计算效率更高的模型。Fernandez Del Rincon等[6]，在此基础上结合用于计算弯曲和根部顺应性的有限元模型和用于计算啮合刚度的接触点近场解析公式，首次将角接触纳入，一定程度上增加了模型分析的准确度。基于计算成本考虑，使用预先计算的啮合刚度替换时变啮合刚度可被接受，但事实证明载荷变化对啮合刚度的影响较大，故均存在计算精度不足的问题。

基于以上研究的不足，最近人们研究了齿廓形状优化对时变啮合刚度和响应特性的影响。Wei等[7]，分别从理论和实验的角度研究了齿轮传动系统的辅助部件对二者的影响，提出了结合齿轮传动系统所有部件动态作用的数值分析模型，该模型的分析结果在额定工况下振动的幅度和频率分布与实验结果吻合。Diez-Ibarbia等[8]，研究了齿廓带齿顶卸载口的齿轮的动态响应特性，其中重点研究了负载分配和摩擦系数的影响。Yu等[9]，分析了齿顶修正对啮合刚度的影响，并指出了在轻载或空载条件下齿轮驱动侧和外侧间啮合刚度的转换关系。

尽管相关研究已经取得了很大的进展，但仍有许多问题尚待解决。其中，很重要的是如混沌运动等非线性接触和传动系统振动及动态响应之间因果关系的模糊。尽管在理论上有许多模型都可以预测这种现象，但实际上，人们在建模的理论上并没有达成共识，从而尚未在计算成本和计算精度之间取得平衡。更重要的是迄今为止，大多数研究人员都忽略了齿轮动态响应的潜在重要组成部分，即单齿的本征振动和非线性接触。即使FEA有所考虑，也常因建模计算成本高忽略它们[10]。

为了解决这些问题，文中开发了一种结合了齿间非线性接触及齿本征振动的高计算效率的简化计算模型b，并通过与常规的全要素模型a对比验证了其正确性。首先模型b通过线性化分析，得到系统的本征振动频率和时变啮合刚度；然后以此为基础，通过非线性仿真检查各种工作状态(如非线性接触和混沌运动)间的过渡；最后得到传动系统的动态响应特性曲线和振动特性曲线。因此该模型涵盖了绝大多数工况，具备了良好的鲁棒性，且在分析成本和精度之间取得了良好的平衡。同时在线性分析和非线性分析中分别考虑了齿的本征振动和非线性接触，较好弥补了前人工作的不足。

1 建模及计算

1.1 齿轮级理论计算建模

对单个齿轮级进行建模，可以表征其旋转动力学特性。建模中，齿轮的惯性特性由每个齿轮的惯性矩(Iti)及等效齿轮惯性矩(Ii)通过扭力弹簧(kti)耦合表示，其中扭力弹簧(kti)对应于每个齿的弯曲、剪切、压缩和齿根顺应性，当指代主动齿轮时，i的索引值为1，当指代从动齿轮时，i的索引值为2。且该弹簧带有扭转阻尼(cti)，它代表与上述参数相关的固有材料阻尼。同时接触齿对通过非线性接触力彼此耦合，用弹簧(kh)表示，齿轮级模型示意图如图1所示。

图1 单齿轮级的简化集总元素模型

在组成传动系统的n个齿轮级模型中，每个齿轮级模型都包括齿轮轮毂的一对自由度(φ1和φ2)，以及每个啮合齿对的另外一对自由度(φt1和φt2)。φi表征每个齿轮毂绕其轴线的旋转，而φt表征每个齿轮齿绕其根部的等效旋转。

模型中涉及的许多几何因素，例如杠杆(Li和Lti)都是可变的，并取决于啮合位置。为了提高计算速度，可以用主动齿轮齿的角位置(φmp)为参考变量，检索并使用预先计算好后储存在查询表中因子的瞬时值。而查询表中，存储了从静态到动态啮合位置的全部几何数据。在总共z个齿的齿轮中，第i齿的φmp=φ1+2πi-z1(从-π转换到π)，其中单齿计算中涉及的几何元素如图2所示。

图2 计算中涉及的几何元素

最后使用C. Spitas[11]提出的数值方法确定接触点，此方法的原理是利用二维啮合方程的一维公式和牛顿-拉夫森迭代方案来定位两个任意齿轮轮廓间的接触点，并且考虑了角接触的情况。因为如果发生明显的齿变形，角接触的作用则不可忽略。通过用恒定半径条件(分别等于主动齿轮或从动齿轮的齿尖半径)替换齿廓之间的相切条件，可以实现对此类情况的估算分析，以提高计算效率。

1.1.1 啮合刚度分析计算

啮合刚度是指齿轮啮合时轮齿抵抗变形的能力。其基本前提是有效区分单个齿的顺应性和源自齿间啮合界面的影响因素。

自由度φti则可作为一个区分量，根据φti得到的模型，能以简化方式合并先前的成分。各成分中，特别需考虑到每个齿的弯曲度(Cb)、剪切力(Cs)、压缩性(Cc)、无量纲反刚度(C)和齿根顺应性(Cf)。对于C的计算公式如下：

C=Ebδ/W

(1)

式中：E为弹性模量；b为齿宽；δ为沿作用线的相对位移；W为法向载荷。

根据ISO 6336-1标准，可使用势能法估算每个单独啮合位置的Cb，Cs和Cc。而对于Cf，则可使用O'Donnel提出的分析解决方案计算[12]。

如上节所述，可将各个量组合在一个等效的扭转弹簧kti中，该弹簧通过可变且同时与啮合位置相关的杠杆Lti加载(如图1)。选择kti是为了使齿的顺应性的单自由度模型解，最接近于所有啮合位置的解析解。对于每个单独的啮合位置，则可使用以下公式计算kti的等效瞬时值：

(2)

将kti与Iti结合使用既可以简化对单齿顺应性的建模，又不会影响加载时对可变啮合刚度特性的表达。计算结果表明，这样既可在加载时表征时变啮合刚度，又可在空载时表征单齿的顺应性。

对于接触顺应性，可通过Weber和Banascheck[13]提出的分析解决方案计算。此方案不仅考虑了接触顺应性对啮合位置的依赖性，还考虑了运行负载，分析计算公式如下：

(3)

式中：ri(i=1,2)为啮合齿在接触点处的局部曲率半径；E、v为材料的弹性特性；δ为压痕；hi(i=1,2)为小齿轮和大齿轮齿的深度。

压痕δ是指沿两点之间作用线的相对位移，其分别位于小齿轮和大齿轮齿的深度hi(i=1,2)中，可以通过φmp查表确定。

对于给定的φmp，可以用空载时的预计算值为基础，并根据空载接触位置的所有角度的微小扰动来计算。将接触位置设置为：

(φ1o,φ2o,φt1o,φt2o)=(φ1,δφ2(φmp)-rφ1,0,0)

(4)

由于φ1和φmp彼此成正比，且两者间的相互作用已包含在预计算中，因此不必考虑φ1的扰动。故δ的瞬时值为：

δ=-{[rφ1+φ2-δφ2(φmp)]L2(φmp)+

φt1Lt1(φmp)+φt2Lt2(φmp)}

(5)

式中：r为传动比，且r=z2/z1。

此外，δφ2表征无负载时的静态传动误差，并被计算为从动齿轮的实际接触位置与标称接触位置的偏差。

1.1.2 惯性特性分析计算

惯性特性可以通过两个不同的转动惯量(Ii和Iti)表示，每个齿轮的转动惯量Ii都表示包括齿轮轮毂和所有齿轮齿的惯性，而每个齿的转动惯量Iti都表示单齿相对于齿根的惯性。计算两个转动惯量时的重叠，使得次级惯性效应(例如离心效应)可以忽略[14]。同时，简化模型合并了与齿本征振动有关的影响因素，并将运动方程以简化形式表示，即：

(6)

1.1.3 啮合接触分析计算

每一个齿都存在内外两侧边缘，故应该对式(6)进行修正，以分别考虑两侧的不同接触状态。而修正方法是对相应项求和，并为每个额外的齿对添加更多方程。此外，由于外侧啮合位置和等效驱动侧啮合位置之间具有轴向对称性，故将外侧接触纳入模型中，且其参数已预先计算。可以通过检查传输符号s的正负进行切换,相关因子s的定义为：

(7)

对于大多数齿轮级而言，因为重叠系数很少超过2，所以仅需考虑已经啮合或最接近啮合位置的三个齿对；同理所考虑的主动齿轮的齿角位置(φmp)也应是啮合齿对的左右相邻齿对预计算出的接触位置与轮心连线的夹角。其中两齿对间接触线对应的变桨点(φ1p)，可以通过比较sφmp和φ1p来判断，即可确定接触位置。综上所述，接触位置可使用以下运动方程组简化计算：

(8)

主动齿轮齿与从动齿轮齿是根据它们相互接触的接近程度来配对，故式(8)可进一步用下式表示：

(9)

基于以上定义可以表征齿轮的啮合齿对随齿轮级的旋转而变化，而这种变化的实质是：

(1) 从一个齿对紧密啮合，到另一对代替它。

(2) 开关量s变化，使啮合齿对的啮合侧(接触位置)发生的转变。

1.1.4 单齿的本征振动及阻尼分析

齿的本征振动表征齿在受到外界持续传动力作用产生的瞬态自由振动，齿在正常或异常状态下都会产生本征振动。齿的本征振动与齿的材料、啮合刚度、惯性特性及形状尺寸密切相关，其中受齿材料特性的影响最大，常用材料的阻尼予以表征。

材料的阻尼往往以两种机制出现。一种阻尼机制通常被模拟为旋转阻尼元件cti与旋转弹簧kti的串联放置，该机制表征单齿本征振动中固有的材料阻尼。其中阻尼比可由下式计算：

(10)

对于钢齿轮可以查表得阻尼比ζ=10-4[15]。

另一种阻尼机制则表征每当一个新齿对进入潜在的啮合区域时，虽然没有开始振动，但仍被赋予振动初始条件的现象。

式(10)反映了一个物理事实，即仅在高速齿轮传动系统中，单齿的本征振动才会被覆盖，而在中低速齿轮系统中则不可忽略。在新模型中，阻尼应该被限制在允许稳态响应发生的范围内，以避免系统阻尼过大，掩盖非线性接触和混沌行为。

1.1.5 线性化分析计算

即使提出的模型为高度非线性化模型，但通过导出等效的线性系统，对于检查系统的本征振动特性并确定与非线性分析相关的许多因素(例如采样率或总模拟时间)都不可或缺[16]。运动方程的线性化是基于以下假设进行的：

(1) 在啮合侧持续接触。

(2) 假定额定载荷恒定，用于计算接触顺应量。

(3) 只考虑了两对齿，其中一对在变桨点接触，而另一对不接触。这两对分别说明了齿轮齿的加载和卸载行为，因此整个齿轮级可由6自由度线性系统表示。

基于这些假设，可以将运动方程式(9)简化为：

(i=1,2)

(11)

(12)

在计算非接触齿对时，将f(δ)=W0= 0且kh=0代入δ并减少等价项后，获得6×6的刚度、阻尼和质量矩阵，并将运动方程简化为线性形式为：

(13)

式中：x=[φ1φt11φt21φt12φt22φ2]T

1.2 传动系统动力学建模

除了完整的齿轮级模型外，还应考虑其他零部件因素的影响，从而保证建成能够预测包含所有对传动系统响应特性有影响因子的分析计算模型。因此，还需对轴和飞轮的扭转动力特性进行了建模。对于轴，可以用单轴的单因素扭转有限元分析来表征，并在其中结合质量、阻尼和弹性特性。就阻尼而言，对于每根轴的本征量，通常选择2%的阻尼比[17]。此值是指一个隔离的轴(即未考虑附加到其上的其他建模元素，如齿轮级和飞轮等)的阻尼比。对于飞轮则使用其惯性矩进行建模即可。

2 分析结果

所建立的集总元素模型，可用于预测图3中描述的典型中等功率传动系统的动力学响应特性。

图3 仿真系统的示意图

该系统的电机被建模为施加在输入轴、齿轮级、输出轴和飞轮上的角速度源，并在其上施加了反向扭矩。为了进一步研究单齿本征振动对其响应特性的影响，则提出了原系统的两种变体。在第一个变体中，使用具有大量小尺寸齿的齿轮(多齿传动系统A)，而在第二个变体中，将齿轮齿数减少而齿的尺寸增大(少齿传动系统B)。但两个变体所选择的齿轮都必须符合下列4个条件：

(1) 两种变体传动系统的传动比相同。

(2) 在两种情况下，等效齿轮的齿轮惯性实际上都相同，且齿端半径在2%的误差范围内相同。

(3) B的输入齿轮的齿数应在17～20之间。

(4) A的模数范围为1.5～2.5 mm，B的模数范围为3.5～5 mm。

这组条件既保证两个传动系统具有明显不同的轮齿动力学特性，又保证了两个传动系统的动力学等效性。两系统齿轮设计参数如表1～4所列。

表1 齿轮参数

表2 齿轮共有的材料属性

表3 飞轮参数

表4 轴的通用参数

2.1 静态分析

对结果的静态分析是为了验证为网格划分一致性而提出的分析公式的正确性。

尽管在参考文献中广泛使用了文中使用的方程式，但几乎没有文献研究类似于少齿传动系统B的，具有低齿数齿轮的齿轮级。在这种情况下，将少齿传动系统B的顺应性分析结果与从有限元模型获得的结果进行比较，以验证对其进行理论分析公式的正确性。

为了直接与分析结果进行比较，有限元模型是基于平面应变建立的二维模型，且假设系统仅有齿对之间发生接触，并使用节点到边缘检测算法来模拟齿接触。同时检查了整个啮合位置范围内施加在齿轮上的支撑和负载，其中包括在主动齿轮内表面上的固定支撑以及在从动齿轮的内表面上的圆柱支撑和扭矩。图4展示了典型啮合位置的几何形状、网格和支撑。并分别对接触齿对(左下)和接触区域(右下)的网格进行了细化。

图4 传动系统B的有限元几何形状和网格

有限元分析的结果可通过下式提取等效值：

(14)

式中：ut为输出齿轮啮合表面上点的切向变形。

由于接触顺应性存在非线性行为，所以总顺应性随施加负载变化而变化，该变化遵循对数趋势，如图5所示。由图5可得，有限元分析与理论分析结果之间的偏差有限，对于1 Nm曲线中最右边的点，最大值为8.3%。除这些点外，所有三个曲线均遵循相似的偏差模式，由有限元结果可预测到，在靠近变桨点的啮合位置处的顺应性值略微升高，而在靠近啮合区域的边缘处的顺应性值略微降低。此外，在提高输出转矩时，两种方法都收敛且得到更相近的结果。这表明，产生这种差异的潜在原因是接触刚度公式在两种方法中不同。在非线性有限元分析中可得，接触刚度受接触区域网格密度的影响。总体而言，理论计算分析的结果与有限元分析结果之间的偏差较小，这表明即使对于低载少齿传动系统B，使用理论分析公式进行建模分析也合理。

图5 三种输出扭矩下动力总成B的理论分析和有限元分析结果对比图

2.2 线性动态分析

动态分析是围绕系统的本征行为进行。通过将描述轴的线性矩阵与描述系统行为的线性矩阵相结合，可以创建等效的刚度和质量矩阵。所得的阻尼本征值和未阻尼本征向量以图形形式表示。图6是当施加载荷为100 Nm(非线性分析中检查的最大负载)时传动系统A的图形形式表示。

图6 传动系统A的线性化系统的无阻尼特征向量矩阵图

图6中，每条线的长度代表等效的元素大小，而深浅不同的线颜色代表相反的元素符号。线长通过平方根函数进行了归一化处理，因此可以表征影响较小的因素，且本征振动频率值对应于系统阻尼比。

在传动系统A的线性化系统无阻尼特征向量矩阵的基础上，使用式(15)进行减少自由度优化，并再次进行分析后，得到动力传动系统B的线性化系统的无阻尼特征向量矩阵如图7所示。

图7 传动系统B的线性化系统的无阻尼特征向量矩阵图

(15)

2.3 非线性动力学仿真分析

非线性分析是通过使用二次开发的Matlab和Simulink代码，并模拟各种情况下的传输响应来进行。其过程首先是使用带有Matlab函数的功能块根据理论分析公式进行非线性运算，例如求解式(3)。然后分别将运动方程式(9)和(11)进行重构，并在Simulink中实现，最终用ode45求解器集成。整个分析检查了包括线性分布在600～6 000 r/min之间的10个随机输入转速值和对数分布在1～100 Nm之间的10个随机输出扭矩值组合得到的100种不同的情况。在两种模型中检查了所有情况，一种使用了完整的非线性运动方程组(9)(模型a)，另一种使用简化的非线性运动方程组(11)(模型b)。在简化模型b中使用了线性化分析得到的时变啮合刚度值，并简化分析了齿的惯性特性。同时其根据静力分析和φ2的值从静态平衡条件得出齿角值，最终使齿角值成为因变量而非自由度。该简化模型b普遍适用于确定齿惯性建模对模拟响应的影响。

2.3.1 仿真分析参数

利用线性动态分析可以估算用于非线性分析的数值参数，如采样率、采样时间和总模拟时间等。

首先根据Nyquist-Shannon定理，并考虑到基于线性分析，则最大本征振动频率为200 kHz，采样率为1000 kHz。

再根据线性化系统特征频率值，选择数据输出的采样时间。最低本征振动频率为500 Hz，最低激励频率是最低转速时的啮合频率fmes·hmin=140 Hz。由于获得了10 Hz的分辨率，因此采样时间取0.1 s就足够用于光谱分析了。

最后确定总的仿真时间，以保证所有瞬态响应现象在输出数据导出之前都完全发生，并被记录。尽管已对初始条件进行了校准，以使齿对处于接触状态，并且所有组件的角速度都已同步，但即使稳态响应是周期性的，也没有确定标准来确定瞬态响应的结束。因此，必须增加间隔时间以保证能将输出数据中观察到的所有非周期性现象都表征出来。根据通过线性分析计算出的稳定时间，设定总时间为所有特征值的最大稳定时间的10倍，即10 s。

2.3.2 仿真分析结果

考虑到在大量不同场景的普适性，则主要通过使用彩色图形来检查仿真结果，以便更容易识别各区域中的不同响应和特征。该分析的主要目标是根据常见响应特征，将响应分为多个类别，并进行分布分析。因为系统长期的接触磨损与混沌响应的出现息息相关，所以轮副的齿之间存在的不间断接触成为了检测最频繁的特征项。将工作区域划分为恒定接触区和间歇接触区，仅能粗略地表示传动系统的工作状态；为了进一步了解传动系统的动力学特性，则需检查更多有代表性的响应，这时便通过使用相位图表述来解决。相位图是根据输出齿轮的旋转相对角度定义的，即te=rφ1+φ2，其本质上描述了主动齿轮和从动齿轮齿的相对位置。本次仿真总共确定了6种不同的响应类型，且每种响应类型表现出的特征及其分布如下所述。

(1) 类型1：连续接触，呈边界清晰的椭圆形凸瓣状，振动不明显，接触损耗小。

类型1的响应特征是齿轮之间仅存在连续接触，且在相位图中收敛得到一个具有两个不同区域的极限循环，而这两个区域分别表示单齿和双齿接触，如图8和9所示。图8的运行条件为21.5 Nm,600 r/min。图9的运行条件 59.9 Nm,600 r/min。

图8 传动系统A的a型1类响应特性

图9 传动系统B的b型1类响应特性

由相位图上可观察可得：两个区域通常都是具有清晰边界的椭圆形凸瓣状，同时还存在相当细微的形状细节。由于较低的啮合刚度允许较高的齿变形值，因此会产生传动误差，所以用于单齿接触的区域应位于等效的双齿接触区域的右侧。对于此类响应，其极端情况则常出现在单齿接触区域的左边界和双齿接触区域的右边界。由接触振动图观察可得：振动模式较为有序且微弱。

(2) 类型2：连续接触，边界存在明显干扰的成花瓣状，振动不明显，接触损耗小。

类型2的是类型1的变体，主要区别是单齿和双齿接触凸角的相对位置，如图10和11所示。图10的运行条件 21.5 Nm,1200 r/min。图11的运行条件为： 59.9 Nm,1800 r/min。

图10 传动系统A的a型2类响应特性

图11 传动系统B的b型2类响应特性

在类型2中，凸角靠得更近，它们的边界之间存在明显干扰，以至于在传动误差值与单齿或双齿接触之间不存在一对一的等价关系。但是单齿接触凸角的位置比双齿接触凸角更靠右。其余与类型1保持一致。

(3) 类型3：以连续接触为主，夹杂间歇接触，振动不明显，接触损耗较小。

类型3的响应是类型2的极端情况，其相位图呈交叉的双环状，且边界处存在明显干扰。类型3在B中出现几率更高且更具代表性，因此针对性以传动系统B的相位图予以表征。如图12所示。图12的运行条件 32 Nm,2800 r/min。

图12 传动系统B的a型3类响应特性

就连续接触间隔的数量和每个单独间隔的持续时间而言，接触损失的程度仍然受到限制。从质量上讲，与类型2相比，具有接触损耗区域的传动误差值范围更小的特点。就单齿和双齿接触凸角而言，它们的相对位置更加偏移，在极端情况下，两个凸角均在其左端接触，出现了失去接触的情况。就振动而言，相比于响应类型1和2更小。

(4)类型4 以间歇接触为主，夹杂连续接触，振动明显，接触损耗大。

因为大量的接触损耗会导致振动模式和整体行为的混乱，所以此类型的响应很少出现。而在本研究中，这种响应类型也仅在动力总成B输入转速为2 400 r/min左右时出现，如图13所示。图13的运行条件为12.9 Nm,2400 r/min。

图13 传动系统B的b型4类响应特性

该类响应的特点为：虽振动模式明显，但在较大的区域里啮合方式为间歇接触啮合。在仔细分析振动模式时发现，只要激励频率不变，则该振动模式则不会随载荷而变化；又因为齿的本征振动频率也不会随负载而变化，所以可推得响应类型4的振动频率与齿的本征振动频率存在谐振。受接触损耗和恢复时间的影响，在加载和未加载的系统中，各齿的本征振动频率之间都存在特定的耦合关系。

(5) 类型5：间歇性接触，没有两侧接触，振动不明显，接触损耗大。

类型5是最常见的喀哒声状反应类型，通常表现为接触力的明显减少和混沌反应的充分发展，如下图14和15所示。图14的运行条件为 4.64 Nm，1 800 r/min；且为了表达清楚，仅描绘了有限数量的啮合循环。图15的运行条件为 7.74 Nm，5 400 r/min；且为了表达清楚，仅描绘了有限数量的啮合循环。

图14 传动系统A的a型5类响应特性

图15 传动系统B的b型5类响应特性

接触损耗的特性是大多数接触损耗间隔都很短。通常，因为两个齿轮的这种不规则的解耦和重新耦合阻止了有序振动模式的发展，传动误差值的区域比发生接触的区域小得多，所以非接触间隔的这种不规则分布表征了混沌响应。

(6) 类型6：间歇性接触，有两侧接触，接触损耗大，振动模式不明显，

类型6是类型5的变体，相比于类型5，其只是发生了双侧接触现象，如图16所示。图16的运行条件为 12.9 Nm,2 400 r/min。

图16 传动系统A的a型6类响应特性

外侧接触出现表明，系统的振动幅度增大。由于在传动系统B中，待啮合的齿之间存在较大的间隙，所以只能在齿数大到存在三个齿对同时接触时，两侧侧接触这种持续时间不会超过总接触时间5%的现象才可能出现，因此通常认为两侧接触仅出现在传动系统A中。

综上所述，上述分析已经以一种特定的相位定义了6种常见响应类型，以便宽泛对应于从有序的连续接触响应到混乱的间歇接触响应的逐渐过渡。通过这些响应类型在操作方面的分布，可以更好地看到这种过渡，如图17、18所示。

图17 传动系统A在理论计算模型的响应类型分布

图18 传动系统B在理论计算模型的响应类型分布

针对传动系统A，在低速重载时，类型1和2相连且类型3广泛的分布在其周围，组成了一个较大的连续接触区域1；在高速或低载时，类型5和6连续且相互穿插的分布在剩余区域形成间歇性接触区域2；整个图中没有出现响应类型4；且非线性化建模确定了原本模糊的区域1和2之间的边界。同时结合图8～17进行分析可推得：一方面，在非超高速时，连续接触仍为主要的接触方式，接触损耗极小。另一方面，由于系统本征振动频率对外部负载没有依赖性，所以输入转速和边界载荷之间没有一一对应的关系。

和传动系统A相比，传动系统B的特征是：由于啮合参数的改变，例如重叠系数的降低，导致间歇接触增多且接触损耗增大，虽然类型3响应数量少量增加，补偿了连续接触区域1的减小，但随着类型1和2的大量减少，连续接触区域1的尺寸总体上还是趋于减小。此外，传动系统B在2 400 r/min的负载范围内产生了罕见的响应类型4。这些间歇接触和高接触损耗响应类型的大量出现表明，传输的临界频率不仅与混乱的响应相关，而且与单齿的本征振动频率相关。

比较传动系统A和传动系统B的分析结果可推得：在这两种情况下，如果不考虑齿的惯性，连续接触区域1的面积都会因响应类型3的少量增加而增大。但是在传动系统A中的增加量微乎其微，而在传动系统B中的增加量则比较明显。该结果表明：齿轮级中齿轮的单齿尺寸和惯性越大，则本征振动频率和其他相关振动对齿轮级的接触响应特性和能耗影响越大。

2.4 计算模型对比验证

如图17、18所示，以多齿传动系统A和少齿传动系统B为载体，分别使用全要素理论计算模型a和简化计算模型b进行分析计算，经过对比可得：两个计算模型响应类型分布结果相近，差异主要零星的出现在特定的高速重载区域；而全要素计算模型a是引用参考文献[7]的数值计算模型，并增加非线性仿真得到，其正确性既已经被前人验证，又在本文2.1中进行了再次验证，因此可以间接证明简化计算模型b的正确性。

3 结 论

文中创新性的建立了一种对齿轮传动系统动态响应及接触振动特性进行分析的简化计算模型。新模型首先通过对齿对啮合的接触特性、啮合刚度、接触位置、阻尼和惯性特性进行综合性的线性化和非线性化优化分析，得到了啮合齿对间的相位和接触特性，再根据啮合齿对间的相位和接触特性，对响应类型进行了分类，并对各响应类型的特点和分布进行了详尽的阐述。其中所有响应类型都充分考虑到了单齿本征振动的影响。

新模型创新性的研究了具有低齿数的齿轮传动系统B，并与常规多齿数传动系统A进行了对比。此外，还在使用了完整运动方程式的全要素计算模型a的基础上简化得到了使用不包含齿的惯性和动力学方程式的简化计算模型b，并最终通过对a和b的线性和非线性建模分析对比，证明了简化计算模型b的正确性。并且在分析过程中证明，简化计算模型b的计算速度相比于全要素计算模型a提升了6.8%，即计算成本下降6.8%。

前人的线性化仿真模型并未充分考虑非线性情况。而新模型包含非线性仿真，并显示出了丰富的非线性响应特性，即连续接触区域和间歇接触区域间的边界几何形状。且其因为受单齿本征振动的影响而变得极其复杂，致使传动系统在输入转速和边界载荷之间没有一一对应的关系。