CN104502096A - 齿轮动态传动误差计算方法 - Google Patents

Abstract

本发明公开了一种齿轮动态传动误差计算方法，其通过由角度编码仪测得的输入轴角速度、由磁粉制动器提供的输出轴阻力矩、齿侧间隙、输入轴和输出轴的角位移等参数值计算出碰撞时间间隔，将所得到的相对位移与齿侧间隙进行比较，直至相对位移的值在适当的区间内并且符合要求，所得到的结果即为精度较高的传动误差。本方案引入了齿侧间隙对传动误差的影响，保证了计算结果的准确性，同时使计算精度和计算时间得到了良好的平衡，适用于对实际齿轮组的测量计算和仿真模拟的计算。

Description

齿轮动态传动误差计算方法

技术领域

本发明涉及机械制造和控制领域，尤其是涉及一种考虑齿侧间隙的齿轮动态传动误差计算方法。

背景技术

在所有的机械传动中，齿轮传动有着传动平稳、传动比精确、效率高、寿命长等特点，使得齿轮传动成为最重要而且应用最广的机械传动。高速，平稳，精密齿轮传动是国家重点科技攻关项目。传动误差已经被人们普遍的接受和认为是齿轮系统振动和噪声的激励源，因为传动误差影响着齿轮传动中的力和速度的变化。同时，齿轮的动态传动误差也是齿轮检测中十分重要的一项内容，精确测试和计算传动误差对齿轮制造和齿轮动力学研究具有重大意义。

中华人民共和国国家知识产权局于2008年10月01公开了公告号为CN101275881的专利文献，名称是小模数齿轮传动误差测量方法，其先将被测齿轮、带动器、第一测角编码器、第一联轴器和第一电机依次连接，再将测量齿轮、柔性联轴器、第二测角编码器、第二联轴器和第二电机依次连接。被测齿轮和测量齿轮单面啮合，并分别由两个电机进行驱动。检测出测量齿轮相对于第二测角编码器的角位移φ2、测量齿轮轴系角位移φ1和被测齿轮轴系角位移φ′，经过计算得到小模数齿轮的传动误差。但是，由于齿轮的加工误差，安装误差，润滑和修行等因素使得齿廓表面相对于理想齿廓位置的偏移会引起齿侧间隙(图1)，导致齿轮在进入和退出啮合时，啮入和啮出点的位置会偏离理论啮合点，产生线外啮合，使啮合齿面间产生碰撞冲击。这种碰撞冲击对高速齿轮传动的平稳性和精密性起着关键的作用。对于精密传动而言，正反转的时候会有微量空行程，而且有冲击，影响传动精度。显然，齿侧间隙直接影响着齿轮系统得传动误差。然而，此方案未考虑齿侧间隙对传动误差的影响，会使最终的计算结果与实际值发生较大偏差。

同时，对于齿轮系统传动误差的数值计算成为了一个日益突出的问题。目前，对于传动误差的计算方法，多停留在应用现有的传统方法，如：Newmark积分法，Runge-Kutta积分法，Gear法和改进的Gear法(Wstiff法和Dstiff法)。上述计算方法对于传动误差的计算方法的侧重点各有不同，有些方法更侧重于对计算精度的把握；也有一些方法拥有较快的计算速度。传统的计算方法通常无法完美的控制好计算时间与计算精度的分配：在精度要求较高的情况下往往带来庞大的计算量，导致计算时间的漫长；在保证计算时间的条件下，又会使计算精度下降以致于无法正确的描述及预测齿轮系统动传动误差。

发明内容

本发明主要是解决现有技术在实际测量过程中因齿侧间隙导致传动误差测量不准确的技术问题，及在理论模拟计算齿轮动态传动误差难以平衡计算精度和计算时间的问题。在考虑到齿侧间隙对于传动误差到影响到前提下，提出了一种将实际测量和理论模拟计算相结合，并能在保证一定的计算精度情况下快速计算齿轮动态传动误差的计算方法。

本发明针对上述技术问题主要是通过下述技术方案得以解决的：一种齿轮动态传动误差计算方法，包括以下步骤：

A、获取参数值，参数值包括由角度编码仪测得的输入轴角速度和由磁粉制动器提供的输出轴阻力矩；

B、计算碰撞时间间隔Δt_i；

C、比较碰撞时间间隔Δt_i和时间步长Δt的大小，如果Δt＜Δt_i，则进入步骤D；如果Δt≥Δt_i，则进入步骤G；

D、以时间步长Δt为步长来计算相对位移x(n+1)，如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于齿侧间隙L，则进入步骤J；如果相对位移x(n+1)的绝对值|x(n+1)|大于齿侧间隙L，则缩小时间步长Δt然后进入步骤E；

E、比较时间步长Δt和最小时间步长t_min，如果时间步长Δt大于或等于最小时间步长t_min，则进入步骤D；如果时间步长Δt小于最小时间步长t_min，则进入步骤F；

F、以最小时间步长t_min为步长计算相对位移 然后进入步骤J；

G、以碰撞时间间隔Δt_i为步长来计算相对位移x(n+1)，如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于L-ε，则进入步骤J；如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于L+ε且相对位移x(n+1)大于L-ε，则进入步骤I；如果相对位移x(n+1)大于L+ε，则减小碰撞时间间隔Δt_i，然后进入步骤H；

H、比较碰撞时间间隔Δt_i和最小时间步长t_min的大小，如果碰撞时间间隔Δt_i大于或等于最小时间步长，则返回步骤G；如果碰撞时间间隔Δt_i小于最小时间步长t_min，则取最小时间步长为步长，即t(n+1)＝t(n)+t_min，然后进入步骤I；

I、计算相对位移使步长数n增大1，再取碰撞时间常数∈_Δt为步长，即t(n+1)＝t(n)+∈_Δt，计算相对速度为： 计算相对位移为：x(n+1)＝x(n)，然后进入步骤J；

J、判断t(n+1)是否大于M*T_f，如果是，则结束计算，如果否，则步长数n增大1然后返回步骤B；

δ为误差小量，M为周期数，T_f为激励周期，t(n)为第n个时间采样点，x(n)为第n个相对位移，符号上方带原点表示对时间的求导，一个圆点表示一阶导数，两个圆点表示二阶导数，sign为符号函数。

最终得到的M个周期的相对位移x(n+1)、相对速度和时间t(n+1)就能良好地反映齿轮的动态传动误差。通过所获得的传动误差，可以对齿轮箱的精度进行评价。

在的对实际齿轮组进行测量时，各项参数值由常规方式测量得到或由用户设定；仿真模拟计算时，参数值由手动设定或者计算机自动设定。参数值包括L、ε、T_f、R_bg、R_bp、θ_g、θ_p、外载荷力矩的均值、从动轮转动惯量、为外载荷力矩的波动部分、ω、t₀、由于油膜引起的在齿轮副接触正面和侧面的啮合刚度、由于油膜引起的在齿轮副接触正面和侧面的啮合阻尼等。

作为优选，步骤B中，所述碰撞时间间隔计算方式如下：

${\mathrm{\Δt}}_i=\frac{-b\±\sqrt{b^2-4\mathrm{ac}}}{2a},$

式中， $a=(R_{\mathrm{bg}}{T}_m+R_{\mathrm{bg}}{T}_p\cos \left(\mathrm{\ω}{t}_0\right)-R_{\mathrm{bg}}\mathrm{kx}\left(t_0\right)-R_{\mathrm{bg}}c\stackrel{\·}{x}\left(t_0\right))/2,b=\stackrel{\·}{x}\left(t_0\right),$ $c=x\left(t_0\right)-\mathrm{Lsign}\left(\stackrel{\·}{x}\left(t_0\right)\right).$ Δt_i取两个值中最小的正实数，R_bg为从动轮的基圆半径，R_bp为主动轮的基圆半径，T_m为外载荷力矩的均值与从动轮转动惯量的比值，T_p为外载荷力矩的波动部分与从动轮转动惯量的比值，ω为外载荷力矩的频率，t₀为系统起始时间，k为由于油膜引起的在齿轮副接触正面和侧面的啮合刚度的和与从动轮转动惯量的比值，c为由于油膜引起的在齿轮副接触正面和侧面的啮合阻尼的和与从动轮转动惯量的比值。

作为优选，计算碰撞事件间隔Δt_i时，如果Δt_i无正实数，则直接判定Δt_i＜Δt。

作为优选，步骤D中，缩小时间步长的具体操作为：将Δt减小Δt²。

作为优选，步骤D中，减小碰撞时间间隔的具体操作为：将Δt_i减小Δt_i ²。

本发明带来的实质性效果是，考虑了齿侧间隙对传动误差的影响，在保证一定的计算精度要求之下，尽量减少计算时间，实现计算精度与计算时间的一种更加优化的分配。

附图说明

图1是本发明的一种齿轮啮合区放大图；

图2是本发明的一种参数测量系统示意图；

图3是本发明的一种主动轮和从动轮按标准中心距装配示意图；

图4是本发明的一种主动轮和从动轮粘着连动时的判断情况的示意图；

图5是本发明的一种发生齿与齿之间碰撞时Δt≥Δt_i的情况的示意图；

图6是本发明的一种利用ε对轮齿之间的碰撞进行再一次精确的情况的示意图；

图7是本发明的一种计算方法流程图；

图中：1电机；2小齿轮；3输入轴；4角度编码仪；5大齿轮；6输出轴；7磁粉制动器。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：图2是对实际齿轮组进行测量的系统，包括电机1、主动轮2、输入轴3、角度编码仪4、从动轮5、输出轴6和磁粉制动器7。主动轮和从动轮按标准中心距装配如图3所示，图1是齿轮啮合区放大图。

齿轮的传动误差(相对位移)为与传动特性有关的沿啮合线方向度量被动轮上的齿廓在实际啮合时所处位置同理想条件下应处位置之间的差值，可以定义为：

x(t)＝R_bpθ_p(t)-R_bgθ_g(t)    (1)

式中，R_bg为从动轮的基圆半径，R_bp为主动轮的基圆半径；θ_g为从动轮的角位移，θ_p为主动轮的角位移，角位移也可以通过角度编码仪进行测量。

角度编码仪：测量图2中的输入轴的角速度

磁粉制动器：给输出轴提供一个阻力矩，该阻力矩为：

${\stackrel{\‾}{T}}_g={\stackrel{\‾}{T}}_m+{\stackrel{\‾}{T}}_p\cos \left(\mathrm{\ωt}\right)---\left(2\right)$

式中，和分别为外载荷力矩的均值和波动部分的力矩；ω为该阻力矩的频率。

当|x(t)|<L时，根据达朗贝尔原理建立该齿轮系统动力学方程：

$I_g{\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_g\left(t\right)-(c_1+c_2)\stackrel{\&CenterDot;}{x}\left(t\right)-(k_1+k_2)x\left(t\right)=-{\stackrel{\&OverBar;}{T}}_g---\left(3\right)$

其中，I_g为大齿轮的转动惯量；k₁和c₁分别为由于油膜引起的在齿轮副接触正面的啮合刚度和啮合阻尼；k₂和c₂分别为由于油膜引起的在齿轮副接触侧面的啮合刚度和啮合阻尼；和分别为外载荷力矩的均值和波动部分的力矩。将公式(3)两边同时除以转动惯量I_g得：

${\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_g\left(t\right)-c\stackrel{\&CenterDot;}{x}\left(t\right)-\mathrm{kx}\left(t\right)=-T_m-T_p\cos \left(\mathrm{\&omega;t}\right)---\left(4\right)$

这里，c＝(c₁+c₂)/I_g,k＝(k₁+k₂)/I_g，和

根据公式(4)，可以得到：

${\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}}_g\left(t\right)=c\stackrel{\&CenterDot;}{x}\left(t\right)+\mathrm{kx}\left(t\right)-T_m-T_p\cos \left(\mathrm{\&omega;t}\right)---\left(5\right)$

从t₀时刻到t时刻，对公式(5)进行积分，可以得到相应的大齿轮的角速度和角位移

$\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_g={\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_g\left(t_0\right)-T_m(t-t_0)-T_p\frac{\sin \left(\mathrm{\&omega;t}\right)-\sin \left(\mathrm{\&omega;}{t}_0\right)}{\mathrm{\&omega;}}+c(x\left(t\right)-x\left(t_0\right))\\ +k(x\left(t_0\right)(t-t_0)+\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{(t-t_0)}^2}{2})\end{array}---\left(6\right)$

$\begin{array}{c}{\mathrm{\&theta;}}_g\left(t\right)={\mathrm{\&theta;}}_g\left(t_0\right)+{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_g\left(t_0\right)(t-t_0)-T_m\frac{{(t-t_0)}^2}{2}+T_p\frac{\cos \left(\mathrm{\&omega;t}\right)-\cos \left(\mathrm{\&omega;}{t}_0\right)}{{\mathrm{\&omega;}}^2}\\ +T_p\frac{\sin \left(\mathrm{\&omega;}{t}_0\right)}{\mathrm{\&omega;}}(t-t_0)+\mathrm{kx}\left(t_0\right)\frac{{(t-t_0)}^2}{2}\\ +k\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{(t-t_0)}^3}{6}+c\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{(t-t_0)}^2}{2}\end{array}---\left(7\right)$

${\mathrm{\&theta;}}_p\left(t\right)={\mathrm{\&theta;}}_p\left(t_0\right)+{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}_p\left(t_0\right)(t-t_0)---\left(8\right)$

将公式(7)与公式(8)带入公式(1)，可得：

$x\left(t\right)=x\left(t_0\right)+\stackrel{\&CenterDot;}{x}\left(t_0\right)(t-t_0)-R_{\mathrm{bg}}\left(\begin{array}{c}{T}_p\frac{\cos \left(\mathrm{\&omega;t}\right)-\cos \left(\mathrm{\&omega;}{t}_0\right)}{{\mathrm{\&omega;}}^2}+T_p\frac{\sin \left(\mathrm{\&omega;}{t}_0\right)}{\mathrm{\&omega;}}(t-t_0)\\ +(\mathrm{kx}\left(t_0\right)+c\stackrel{\&CenterDot;}{x}\left(t_0\right)-T_m)\frac{{(t-t_0)}^2}{2}+k\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{(t-t_0)}^3}{6}\end{array}\right)---\left(9\right)$

$\stackrel{\&CenterDot;}{x}\left(t\right)=\stackrel{\&CenterDot;}{x}\left(t_0\right)-R_{\mathrm{bg}}\left(\begin{array}{c}c(x\left(t\right)-x\left(t_0\right))+k(x\left(t_0\right)(t-t_0)+\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{(t-t_0)}^2}{2})\\ -T_m(t-t_0)-T_p\frac{\sin \left(\mathrm{\&omega;t}\right)-\sin \left(\mathrm{\&omega;}{t}_0\right)}{\mathrm{\&omega;}}\end{array}\right)---\left(10\right)$

当|x(t)|＝L时，齿轮系统发生弹性碰撞，此时系统动力学方程为：

$x(t+{\&Element;}_{\mathrm{\&Delta;t}})=x\left(t\right),\stackrel{\&CenterDot;}{x}(t+{\&Element;}_{\mathrm{\&Delta;t}})=-e\stackrel{\&CenterDot;}{x}\left(t\right)---\left(11\right)$

其中，e为碰撞系数，碰撞时间常数∈_Δt为非常小的一个正数。在碰撞时，传动误差可以表示为：

$x(t_0+{\mathrm{\&Delta;t}}_i)=\mathrm{L\; sign}\left(\stackrel{\&CenterDot;}{x}\left(t_0\right)\right)---\left(12\right)$

这里，Δt_i＝t_i-t₀为碰撞时间间隔；若则表示碰撞发生在轮齿的正面，此时，从动轮在主动轮的推动下转动。若则表示碰撞发生在轮齿的侧面，此时，主动轮在从动轮的推动下转动。

将公式(9)结合公式(12)得：

$\begin{array}{c}\mathrm{L\; sign}\left(\stackrel{\&CenterDot;}{x}\left(t_0\right)\right)=x\left(t_0\right)+\stackrel{\&CenterDot;}{x}\left(t_0\right){\mathrm{\&Delta;t}}_i\\ -R_{\mathrm{bg}}\left(\begin{array}{c}{T}_p\frac{\cos \left(\mathrm{\&omega;t}\right)-\cos \left(\mathrm{\&omega;}{t}_0\right)}{{\mathrm{\&omega;}}^2}+T_p\frac{\sin \left(\mathrm{\&omega;}{t}_0\right)}{\mathrm{\&omega;}}{\mathrm{\&Delta;t}}_i\\ +(\mathrm{kx}\left(t_0\right)+c\stackrel{\&CenterDot;}{x}\left(t_0\right)-T_m)\frac{{\mathrm{\&Delta;t}}_i^2}{2}+k\stackrel{\&CenterDot;}{x}\left(t_0\right)\frac{{\mathrm{\&Delta;t}}_i^3}{6}\end{array}\right)\end{array}---\left(13\right)$

对公式(13)左右两边同时在小量Δt_i附件进行2阶泰勒公式展开，并忽略比高阶的无穷小量，得：

$\mathrm{L\; sign}\left(\stackrel{\&CenterDot;}{x}\left(t_0\right)\right)\&ap;x\left(t_0\right)+\stackrel{\&CenterDot;}{x}\left(t_0\right){\mathrm{\&Delta;t}}_i-\frac{R_{\mathrm{bg}}}{2}(\mathrm{kx}\left(t_0\right)+c\stackrel{\&CenterDot;}{x}\left(t_0\right)-T_m-T_p\cos \left(\mathrm{\&omega;}{t}_0\right)){\mathrm{\&Delta;t}}_i^2---\left(14\right)$

从上式可以求出：

${\mathrm{\&Delta;t}}_i=\frac{-b\&PlusMinus;\sqrt{b^2-4\mathrm{ac}}}{2a}---\left(15\right)$

式中， $a=(R_{\mathrm{bg}}{T}_m+R_{\mathrm{bg}}{T}_p\cos \left(\mathrm{\&omega;}{t}_0\right)-R_{\mathrm{bg}}\mathrm{kx}\left(t_0\right)-R_{\mathrm{bg}}c\stackrel{\&CenterDot;}{x}\left(t_0\right))/2,b=\stackrel{\&CenterDot;}{x}\left(t_0\right),$ $c=x\left(t_0\right)-\mathrm{Lsign}\left(\stackrel{\&CenterDot;}{x}\left(t_0\right)\right).$ 这里根据物理意义，Δt_i只取其中最小的正实数。若Δt_i无最小正实数，则认为Δt<Δt_i。

如图7所示，将Δt_i与Δt进行比较来确定下一个时间步长为多少，如果是Δt<Δt_i，则以时间步长Δt来计算相对位移x(n+1)。若|x(n+1)|<L，则计算下一时刻的时间Δt_i。若|x(n+1)|>L，则说明相对位移大于齿侧间隙，此时应该按照公式

△t＝△t-△t²    (16)

缩小时间步长Δt。若公式(16)缩小后的时间步长大于系统所设置的最小时间步长t_min，则重新计算相对位移x(n+1)。反之，系统以最小时间步长t_min为步长计算相对位移x(n+1)。此时系统大小齿轮发生粘动现象(如图4所示)，此时根据其物理意义，可以定义

$x(n+1)=\mathrm{Lsign}\left(\stackrel{\&CenterDot;}{x\left(t\left(n\right)\right)}\right),\stackrel{\&CenterDot;}{x}(n+1)=0---\left(17\right)$

如果是Δt≥Δt_i(如图5所示)，则用Δt_i作为下一个步长去计算下一次的相对位移与速度。同样，这里我们要对x(n+1)进行判断，如果|x(n+1)|>L+ε，则证明这里的相对位移x(n+1)已经超出了精确范围，如图6所示，这里就需要减小Δt_i，来保证精确到允许的误差范围内，减小Δt_i的公式如下：

${\mathrm{\&Delta;t}}_i={\mathrm{\&Delta;t}}_i-{\mathrm{\&Delta;t}}_i^2---\left(18\right)$

以公式(18)中的Δt_i为步长，重新计算x(n+1)与如果此时|x(n+1)|>L+ε，重复公式(18)，直到Δt_i<t_min发生时，系统将以t_min为步长计算相对位移和相对速度。此时，系统发生碰撞，根据碰撞的实际物理意义，可以令：

$x(n+1)=\mathrm{Lsign}\left(\stackrel{\&CenterDot;}{x\left(t\left(n\right)\right)}\right),---\left(19\right)$

为保障系统在碰撞前后，其相对位移值相等，这里引入一个时间小量，即碰撞时间常数∈_Δt，使得系统以∈_Δt为步长，即：

t(n+1)＝t(n)+∈_△t    (20)

重新计算一次相对位移和相对速度，其碰撞后的相对速度和相对位移为：

$\stackrel{\&CenterDot;}{x}(n+1)=-e\stackrel{\&CenterDot;}{x}\left(n\right),x(n+1)=x\left(n\right)---\left(21\right)$

若不满足公式|x(n+1)|>L+ε，且满足|x(n+1)|>L-ε,系统发生碰撞，此时重复公式(19)，(20)和(21)。

若不满足|x(n+1)|>L+ε,也不满足|x(n+1)|>L-ε，说明系统没有碰撞，此时重新计算位移和速度。

这样不断循环，可以记录下大量的x(n+1)，t(n+1)等数据。在其整个运算过程中，若t(n+1)<M*T_f，则把其计算得到的相对位移和相对速度代入公式(15)中重新计算Δt_i。若t(n+1)>M*T_f则终止运算，即可得到M个周期的相对位移(传动误差)，相对速度(传动误差的速度)和时间的数据。M为激励的周期数，T_f为激励周期，激励由磁粉制动器提供。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

尽管本文较多地使用了相对位移、碰撞时间间隔等术语，但并不排除使用其它术语的可能性。使用这些术语仅仅是为了更方便地描述和解释本发明的本质；把它们解释成任何一种附加的限制都是与本发明精神相违背的。

Claims (5)

1.一种齿轮动态传动误差计算方法，其特征在于，包括以下步骤：

A、获取参数值；

B、计算碰撞时间间隔Δt_i；

C、比较碰撞时间间隔Δt_i和时间步长Δt的大小，如果Δt＜Δt_i，则进入步骤D；如果Δt≥Δt_i，则进入步骤G；

D、以时间步长Δt为步长来计算相对位移x(n+1)，如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于齿侧间隙，则进入步骤J；如果相对位移x(n+1)的绝对值|x(n+1)|大于齿侧间隙L，则缩小时间步长Δt然后进入步骤E；

E、比较时间步长Δt和最小时间步长t_min，如果时间步长Δt大于或等于最小时间步长t_min，则进入步骤D；如果时间步长Δt小于最小时间步长t_min，则进入步骤F；

F、以最小时间步长t_min为步长计算相对位移 $\stackrel{.}{x}(n+1)=0,$ 然后进入步骤J；

G、以碰撞时间间隔Δt_i为步长来计算相对位移x(n+1)，如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于L-ε，则进入步骤J；如果相对位移x(n+1)的绝对值|x(n+1)|小于或等于L+ε且相对位移x(n+1)大于L-ε，则进入步骤I；如果相对位移x(n+1)大于L+ε，则减小碰撞时间间隔Δt_i，然后进入步骤H；

H、比较碰撞时间间隔Δt_i和最小时间步长t_min的大小，如果碰撞时间间隔Δt_i大于或等于最小时间步长，则返回步骤G；如果碰撞时间间隔Δt_i小于最小时间步长t_min，则取最小时间步长为步长，即t(n+1)＝t(n)+t_min，然后进入步骤I；

I、计算相对位移使步长数n增大1，再取碰撞时间常数∈_Δt为步长，即t(n+1)＝t(n)+∈_Δt，计算相对速度为： 计算相对位移为：x(n+1)＝x(n)，然后进入步骤J；

J、判断t(n+1)是否大于M×T_f，如果是，则结束计算，如果否，则步长数n增大1然后返回步骤B；

ε为误差小量，M为周期数，T_f为激励周期，t(n)为第n个时间采样点，x(n)为第n个相对位移。

2.根据权利要求1所述的齿轮动态传动误差计算方法，其特征在于，步骤B中，所述碰撞时间间隔计算方式如下：

$\mathrm{\&Delta;}{t}_i=\frac{-b\&PlusMinus;\sqrt{b^2-4\mathrm{ac}}}{2a},$

式中， $a=(R_{\mathrm{bg}}{T}_m+R_{\mathrm{bg}}{T}_p\cos \left({\mathrm{\&omega;t}}_0\right)-R_{\mathrm{bg}}\mathrm{kx}\left(t_0\right)-R_{\mathrm{bg}}c\stackrel{.}{x}\left(t_0\right))/2,$ $b=\stackrel{.}{x}\left(t_0\right),$ $c=x\left(t_0\right)-\mathrm{Lsign}\left(\stackrel{.}{x}\left(t_0\right)\right),$ Δt_i取两个值中最小的正实数，R_bg为从动轮的基圆半径，R_bp为主动轮的基圆半径，T_m为外载荷力矩的均值与从动轮转动惯量的比值，T_p为外载荷力矩的波动部分与从动轮转动惯量的比值，ω为外载荷力矩的频率，t₀为系统初始时间，k为由于油膜引起的在齿轮副接触正面和侧面的啮合刚度的和与从动轮转动惯量的比值，c为由于油膜引起的在齿轮副接触正面和侧面的啮合阻尼的和与从动轮转动惯量的比值。

3.根据权利要求2所述的齿轮动态传动误差计算方法，其特征在于，计算碰撞事件间隔Δt_i时，如果Δt_i无正实数，则直接判定Δt_i＜Δt。

4.根据权利要求1或2或3所述的齿轮动态传动误差计算方法，其特征在于，步骤D中，缩小时间步长的具体操作为：将Δt减小Δt²。

5.根据权利要求1或2或3所述的齿轮动态传动误差计算方法，其特征在于，步骤D中，减小碰撞时间间隔的具体操作为：将Δt_i减小