本发明属于机械动力学技术领域，具体涉及到一种齿轮副啮合特性 分析方法。
背景技术：
齿轮传动中，啮合刚度激励是齿轮啮合的主要动态激励之一，对啮 合特性的准确评估对于直齿轮副动力学性能预测及改进设计至关重要。 目前，现有的直齿轮副啮合特性分析方法具有以下特点：
1)有限元方法
有限元方法在模拟齿轮副接触过程时，通常建立接触单元来模拟真 实的接触关系。这种方法需要足够细密的接触区域网格，才能保证收敛 性。因为涉及大规模的非线性迭代解算，导致该方法的效率很低；但该 方法得出的结果较为准确。
2)解析方法
解析方法中，材料力学方法将轮齿视为变截面悬臂梁，利用变形法 或者势能法进行轮齿变形的求解，在求解轮齿刚度时具有很高的精度， 但是基体刚度的求解精度尚待提高，并且材料力学方法很难考虑真实的 基体结构(如腹板结构、减重孔)对刚度的影响。而经验公式法具有计 算简便的特点，但是其精度往往难以得到保证。整体上，解析方法具有 高效率的特点，但精度不如有限元方法。
另一方面，齿轮传动中，弯曲疲劳是是齿轮的主要失效形式之一。 获得准确的齿根的弯曲应力，是预测齿根弯曲疲劳失效路径及预测齿轮 副服役时间的主要依据。目前，现有的直齿轮副弯曲应力分析方法具有 与前述的啮合特性分析相似的特点。
上述缺陷是本领域技术人员期望克服的。
技术实现要素：
(一)要解决的技术问题
本发明提供一种齿轮副啮合特性分析方法，已解决目前齿轮副啮合 特性分析效率偏低的问题。(二)技术方案
为了达到上述目的，本发明采用的主要技术方案包括：
一种齿轮副啮合特性分析方法，包括以下步骤：
步骤S10：针对获取的齿轮副的本地有限元模型，在啮合点周围区域 刚性化和/或约束主动轮和从动轮的齿根圆附近节点的全部自由度，得到 第一总体刚度矩阵
根据所述第一总体刚度矩阵确定所述齿轮副的轮齿刚度ktooth(τ)；
步骤S20：针对获取的齿轮副的本地有限元模型，对轮齿区域刚性化 和/或约束主动轮和从动轮的内孔附近节点的全部自由度，得到第二总体 刚度矩阵
根据所述第二总体刚度矩阵分别确定主动轮的基体刚度kfp(τ)和从动轮的基体刚度kfg(τ)；
步骤S30：根据所述齿轮副的轮齿刚度ktooth(τ)、主动轮的基体刚度kfp(τ)和从动轮的基体刚度kfg(τ)，确定所述齿轮副的时变啮合刚度k(τ)。
进一步地，所述的分析方法，在步骤S10或步骤S20之前，还包括：
步骤S100：根据获取的直齿齿轮副的几何参数及运行参数，生成建 模文件；
步骤S200：根据所述建模文件，由有限元建模工具软件建立所述直 齿齿轮副的通用有限元模型，并输出描述所述有限元模型的单元-节点文 件；
步骤S300：根据获取的所述有限元模型的单元-节点文件，由解算工 具软件建立所述直齿齿轮副的本地有限元模型。
进一步地，所述的分析方法，所述步骤S10中，约束主动轮和从动 轮的齿根圆附近节点的全部自由度，包括：
约束位于(1±α)*rf区域内的有限元单元的各节点沿x方向和y方向 的直线移动自由度，其中，rf为齿根圆半径，α为0与0.05之间的常数， 所述有限元单元为平面四边形等参数单元。
进一步地，2所述的分析方法，所述步骤S10中，在啮合点周围区域 刚性化，包括：
将以啮合点为圆心、以β*m为半径内的有限元单元的弹性模量的数 值调整为原始弹性模量的γ倍，其中，β为大于0.05、小于0.4的常数， m为齿轮模数；γ为大于500且小于1000的常数。
进一步地，所述的分析方法，所述步骤S20中，约束主动轮和从动 轮的内孔附近节点的全部自由度，包括：
约束位于(1+b)*ri区域内的有限元单元的各节点沿x方向和y方向 的直线移动自由度，其中，ri为主动轮或从动轮的内孔半径，b为0与1 之间的常数，所述有限元单元为平面四边形等参数单元。
进一步地，所述的分析方法，所述步骤S20中，对轮齿区域刚性化， 包括：
将位于齿根圆与齿顶圆之间区域内的有限元单元的弹性模量的数值 调整为原始弹性模量的a倍，其中，a为大于500且小于1000的常数。
进一步地，所述的分析方法，所述步骤S10中，确定所述齿轮副的 轮齿刚度ktooth(τ)，包括：
根据第一公式确定轮齿静态传递误差Er(τ)，所述第一公式为：
其中，kti为参与啮合的第i对轮齿的啮合刚度，i＝1,2,3；Epmax为最大 齿廓偏差；Ep1(τ)和Ep2(τ)分别为参与啮合的两对轮齿的齿廓偏差；Sa为 接近距离；Sr为分离距离；F为总啮合力。
进一步地，所述的分析方法，所述直齿齿轮副的几何参数包括：修 型长度和修形量。
进一步地，所述的分析方法，还包括：
步骤S40：针对获取的齿轮副的本地有限元模型，约束主动轮和从动 轮的内孔附近节点的全部自由度；
根据载荷分配系数在啮合点施加啮合力，进行静力学求解后获得各 个单元的节点位移；
根据所述节点位移，将不同单元的同一节点的应力作平均化处理， 并确定各节点处的von Mises应力σ(τ)。
进一步地，所述的分析方法，还包括：
确定主动轮受拉侧过渡曲线上应力最大值作为齿根弯曲应力σb(τ)0；
根据第二公式确定修正系数λ，所述第二公式为：
其中，τ为无量纲啮合时间；
根据第三公式确定修正后的齿根弯曲应力σb(τ)，所述第三公式为：
(三)有益效果
本发明的有益效果是：本发明提供的齿轮副啮合特性分析方法基于 有限元模型采用解析方法分别分析齿轮副的整体变形及局部接触变形， 并得到时变啮合刚度及齿根弯曲应力。该方法避免了非线性接触迭代， 具有较高的计算效率和较高的计算精度；该方法能考虑直齿轮副的齿顶 修形、延长啮合和复杂齿基。
与现有技术相比，本发明提供的齿轮副啮合特性分析方法的计算时 间远小于ANSYS有限元方法；针对时变啮合刚度的计算误差较小，其精 度可以满足一般情况下对啮合刚度的精度要求；针对弯曲应力的计算误 差较小，其精度可以满足一般情况下对弯曲应力的精度要求。
附图说明
图1为本发明一个实施例的齿轮副啮合特性分析方法的流程图；
图2为一对具有较大修形量的直齿轮副的有限元模型示意图；
图3为一对具有较大修形量的直齿轮副的轮齿变形求解示意图；
图4为一对具有较大修形量的直齿轮副的提前啮合和延迟啮合的示 意图；
图5为本发明实施例的方法分析一对具有较大修形量的直齿轮副的 基体刚度时的示意图；
图6为本发明实施例的方法与ANSYS分析方法分别分析不同扭矩下 一对具有较大修形量的直齿轮副的时变啮合刚度的对比图；
图7为本发明实施例的方法与ANSYS分析方法分别分析一对不修形 的直齿轮副的啮合特性的对比图；
图8为本发明实施例的方法与ANSYS分析方法分别分析一对不修形 的直齿轮副得到的弯曲应力云图。
具体实施方式
为了更好的解释本发明，以便于理解，下面结合附图，通过具体实 施方式，对本发明作详细描述。
为便于理解本发明，对齿轮啮合理论进行如下介绍：根据齿轮啮合 理论，在齿轮副啮合的一个啮合周期T内的任一时刻，最多有两对轮齿 同时参与啮合；在齿轮副啮合的一个啮合周期T内，主动轮或从动轮在 其啮合区，最多有相邻的3对齿轮先后进入和啮合区。也即，存在提前 啮合和延迟啮合的情形。
在计算齿轮副的时变啮合刚度时，需要先分别计算主动轮和从动轮 的轮齿刚度。在齿轮啮合理论中，在单齿啮合区，主动齿轮的轮齿刚度 与从动齿轮的轮齿刚度为串联的关系。在双齿啮合区，主动齿轮的轮齿 刚度与从动齿轮的轮齿刚度为串联的关系，在两对啮合轮齿之间，两对 啮合轮齿的轮齿刚度为并联关系。
具体地，当齿轮副处于双齿啮合区，齿轮副的轮齿刚度如下式所示：
当齿轮副处于单齿啮合区，齿轮副的轮齿刚度如下式所示：
其中，下标的第一个数字的含义为：1代表主动轮，2代表从动轮；
下标的第二个数字的含义为：1代表第一对啮合轮齿，2代表第二对 啮合轮齿；τ为无量纲啮合时间，其上限值为齿轮的重合度，也即无量纲 啮合周期T。
需要说明的是，无量纲啮合周期T的确定方法为本领域技术人员所 公知，这里不再赘述。
另外，需要说明的是，以下内容中，与无量纲啮合时间τ有关的变量 均在其下标中给出了标识τ。
需要说明的是，本发明的齿轮副啮合特性分析方法主要针对直齿轮 啮合情况进行分析。
本发明提供的齿轮副啮合特性分析方法基于有限元模型采用解析方 法分别分析齿轮副的整体变形及局部接触变形，并得到时变啮合刚度及 齿根弯曲应力。该方法避免了非线性接触迭代，具有较高的计算效率和 较高的计算精度；该方法能考虑直齿轮副的齿顶修形、延长啮合和复杂 齿基。
本发明提供的齿轮副啮合特性分析方法的计算时间远小于ANSYS 有限元方法；针对时变啮合刚度的计算误差较小，其精度可以满足一般 情况下对啮合刚度的精度要求；针对弯曲应力的计算误差较小，其精度 可以满足一般情况下对弯曲应力的精度要求。
具体地，本发明实施例的齿轮副啮合特性分析方法，包括如下步骤：
步骤1：获取齿轮副基本参数，利用ANSYS软件参数化编程语言(APDL)在ANSYS中建立齿轮副有限元模型，其中，齿轮副基本参数 包括几何参数和、材料特性参数，在整个分析过程中，均认为为为时不 变常数，具体如表1所列。需要说明的是，在ANSYS中建立的齿轮副有 限元模型中，采用的是二维平面有限元单元，更具体地，所有单元都采 用的平面四边形等参数单元。
步骤2：在ANSYS中导出由节点坐标和单元节点号表示的有限元模 型。具体地，该有限元模型为单元-节点文件，是一个多维矩阵。该多维 矩阵可以采用为本领域技术人员所公知的格式存储在txt文本中。
尽管ANSYS具有组集功能，但步骤2中输出的有限元模型并未进行 过组集操作。
步骤3：在MATLAB科学计算软件中导入该有限元模型，在读取出 节点坐标和单元节点号后，组集总体刚度矩阵，记为K；并在MATLAB 中生成本地有限元模型。在MATLAB中生成的本地有限元模型如图2所 示。
针对在步骤2中并未进行过组集操作的有限元模型进行组集操作， 可以获得初始总体刚度矩阵K。
在MATLAB中生成本地有限元模型，并在MATLAB的显示界面中 显示该本地有限元模型，可以借助利用人工校验，检验导入的有限元模 型是正确的。
需要说明的是，总体刚度矩阵通常为2n*2n的方阵，其中，n为节点 总数目。如，某对啮合齿轮副的有限元模型对应的n大约为3000。
步骤4：确定齿轮副的时变啮合刚度。时变啮合刚度根据轮齿刚度与 齿基刚度来获得。
具体地，齿基刚度包括从动轮的基体刚度和主动轮的基体刚度。
步骤41：确定该齿轮副的轮齿刚度。
第一，约束步骤3中的本地有限元模型中主动轮和从动轮的齿根圆 附近节点的全部自由度。
具体地，在齿根圆附近，约束位于(1±α)*rf区域内的各平面四边形 等参数单元的各节点沿x方向和y方向的直线移动自由度。其中，rf为齿 根圆半径，α为分析系数，为0与0.05之间的常数。
通过向这些节点施加位移约束，可以在一定程度上消除总体刚度矩 阵的奇异性，提高数值方法的分析精度，减小分析误差。
优选地，α为0.01，也即约束位于1.01*rf到0.99*rf区域内的各平面 四边形等参数单元的各节点沿x方向和y方向的直线移动自由度。在一 个示例中，可以将分析误差从20％左右降低至2％。
第二，为了防止力载荷作用于啮合点时，啮合点处产生严重的局部 变形而偏离真实的物理情况，将以啮合点为圆心、以β*m为半径内的全 部有限元单元进行刚性化处理。其中，β为大于0.05、小于0.4的常数， m为齿轮模数。
优选地，β为0.2。
具体地，将这些有限元单元的弹性模量E的数值调整为原始弹性模 量(通常为材料的物理属性)的γ倍，即完成了刚性化。通常，γ为大于 500且小于1000的常数。
优选地，γ为1000，即将这些有限元单元的弹性模量E的数值调整 为原始弹性模量的1000倍。
当集中力施加在啮合点时，会出现失真的局部接触变形。目前常见 处理的方法是建立一个整体有限元模型和一个局部有限元模型。具体地， 在局部有限元模型上施加反向的单位力，得到啮合点的局部变形量，利 用这个局部变形量来抵消整体有限元模型的失真的局部接触变形量，从 而得到准确的整体变形量。
本发明实施例的处理方法是，针对已经建立的整体有限元模型，采 用局部刚性化相应的有限元单元的方式来防止产生局部过大变形；不需 要另外建立局部有限元模型；因此，操作简单、精度已能满足使用要求。
经过约束齿根圆附近节点的自由度以及进行啮合点局部刚性化处理 之后，更新初始总体刚度矩阵K，记这时的总刚度矩阵为
更新总体刚度矩阵的方法如组集等，为本领域技术人员公知，这里 不再赘述。
随着啮合过程推进，啮合点局部刚性化的位置在持续变化，因此， 总刚度矩阵是与无量纲啮合时间τ有关的变量。
需要说明的是，约束齿根圆附近节点的自由度与啮合点局部刚性化 处理这两个步骤没有先后关系，可以根据需要，灵活安排执行的顺序。
第三、获取啮合齿轮副的轮齿刚度
以齿面上一个啮合点处的情形为例，对齿面啮合过程的力与位置的 关系描述如下：
如图3(a)所示，将力载荷F作用于一个啮合点(记这个啮合点在整个 啮合区间内对应的无量纲啮合时间为(τ)从而形成对应于该啮合点处的啮 合压力角β(τ)的第一载荷列阵
随着啮合过程推进，啮合点在齿面的位置持续变化，啮合压力角β(τ)持 续变化，因此，力施加的位置持续变化，也即力施加的节点持续变化， 因此，第一载荷列阵是与无量纲啮合时间τ有关的变量。
求解如式(1)中的线性方程组，得到该啮合齿轮副的第一位移列阵
从第一位移列阵中提取出参与啮合的主动轮的轮齿和从动轮的轮 齿在啮合点处沿啮合线方向的位移upi(τ)和ugi(τ)；
本发明实施例的方法中，设定在啮合点处沿啮合线方向的位移与力 载荷服从线性关系。第一载荷列阵中各元素绝对数值的大小，只影响 该处的齿面位移大小，并不会不影响轮齿刚度。为简单起见，第一载荷 列阵可以采用单位载荷。
进一步地，以第一载荷列阵为单位载荷时确定的主动轮的轮齿和 从动轮的轮齿在啮合点处沿啮合线方向的位移upi(τ)和ugi(τ)，根据式(2)， 得到参与啮合的主动轮的轮齿刚度kpi(τ)和参与啮合的从动轮的轮齿刚度 kgi(τ)：
式(2)中，i表示参与啮合的第i对(i＝1,2,3)轮齿。
作为示例，图3(b)示出了在图3(a)中示出的啮合点处，逐点进 行局部刚性化后主动轮的轮齿位移对比。从图3(b)中可以看出，局部 刚性化有效地抑制了接触点的失真变形。
根据式(3)，在该啮合点处，确定参与啮合的第i对轮齿的啮合刚 度为：
其中，khi为在该啮合点处，参与啮合的第i对轮齿的接触刚度，可以 根据式(4)确定：
式(4)中，E为弹性模量；L为齿宽；F为总啮合力；lsri(τ)为在该 啮合点处，参与啮合的第i对轮齿的载荷分配系数。
其中，总啮合力F＝T1/rb1，rb1为主动轮的基圆半径，T1为施加在主 动轮上的扭矩；或，总啮合力F＝T2/rb2，rb2为从动轮的基圆半径，T2为施加在从动轮上的扭矩。
一对啮合的齿轮副中，主动轮和从动轮的弹性模量通常是相同的。 本申请不涉及主动轮和从动轮的弹性模量不同时的情景。
具体地，若该啮合点位于在单齿啮合区，则这时只有1对轮齿参与 啮合，参与啮合的该对轮齿的载荷分配系数为lsri(τ)为1；若该啮合点位 于在双齿啮合区，则参与啮合的第1对轮齿和第2对轮齿的载荷分配系 数为lsri(τ)均可以选择为0.5。
进一步地，根据式(5)，确定在该啮合点处，轮齿静态传递误差Er(τ)：
其中，kti为参与啮合的第i对轮齿的啮合刚度，即kt1、kt2，kt3；Epmax为最大齿廓偏差；Ep1(τ)和Ep2(τ)分别为参与啮合的两对轮齿的齿廓偏差； Sa为接近距离，两个即将进入啮合的轮齿在啮合线方向上的距离；Sr为 分离距离，两个即将退出啮合的轮齿在啮合线方向上的距离；F为总啮合 力。
在如图4所示，在整个啮合周期内，啮合齿轮副依次为单齿啮合、 双齿啮合和三齿啮合的状态，从而导致较大的冲击振动。具体地，延迟 啮合时，为1、2对齿双齿啮合；提前啮合时，为2、3对齿双齿啮合。
在针对性地对啮合齿轮副修形之后，延迟啮合和提前啮合现象变得 微弱甚至不再发生延迟啮合和提前啮合现象。
具体地，轮齿无修形时，齿廓偏差为零；轮齿修形时，最大齿廓偏 差、齿廓偏差根据修型长度和修形量确定，这里不再赘述。
需要说明的是，工程上并不对单齿啮合区进行修形。
进一步地，根据式(6)，确定啮合齿轮副的轮齿刚度ktooth(τ)：
进一步地，根据式(7)，确定在双齿啮合时参与啮合的两对轮齿间 的啮合力f1(τ)和f2(τ)：
进一步地，根据式(8)，确定双齿啮合时参与啮合的两对轮齿间的载 荷分配系数lsr1(τ)，lsr2(τ)：
步骤42：确定该啮合齿轮副的基体刚度。
第一，将轮齿刚性化。
具体地，对步骤3中的本地有限元模型中的轮齿上全部的有限元单 元(即位于齿根圆与齿顶圆之间的区域内的有限元单元)进行刚性化处 理。具体地，将这些有限元单元的弹性模量E的数值调整为原始弹性模 量(通常为材料的物理属性)的a倍，通常，a为大于500且小于1000 的常数。
优选地，a为1000，也即将这些有限元单元的弹性模量E的数值调 整为原始弹性模量的1000倍。
通过轮齿刚性化步骤，将轮齿变形与基体变形分离开，使得分析得 到的轮齿变形更准确。
在步骤3中确定的总体刚度矩阵K的基础上，记录轮齿刚性化之后 的总刚度矩阵为
第二，约束内孔节点。
分别约束主动轮和从动轮的内孔处节点的全部自由度。
具体地，在内孔处，约束位于(1+b)*ri区域内的各平面四边形等参 数单元的各节点沿x方向和y方向的直线移动自由度。其中，ri为表1中 列出的主动轮或从动轮的内孔半径。
通过向内孔附近的节点施加位移约束，可以在一定程度上消除总体 刚度矩阵的奇异性，提高数值方法的分析精度，减小分析误差。
优选地，b为0.01，也即约束位于1.01*ri区域内的各平面四边形等 参数单元的各节点沿x方向和y方向的直线移动自由度。在一个示例中， 可以将分析误差从20％左右降低至2％。
在步骤3中确定的总体刚度矩阵K的基础上，记录对内孔节点约束 后的总刚度矩阵为
需要说明的是，约束内孔附近节点的自由度与轮齿刚性化处理这两 个步骤没有先后关系，可以根据需要，灵活安排执行的顺序。
经过约束内孔节点的自由度以及进行轮齿刚性化处理之后，再次更 新总体刚度矩阵，记这时的总刚度矩阵为
第三、获取基体刚度。
如图5所示，向齿面施加力载荷，记第二载荷列阵为
求解式(9)中的线性方程组，得到该啮合齿轮副的第二位移列阵
从第二位移列阵中提取出主动轮和从动轮的啮合点沿啮合线方 向的位移ufp(动轮和ufg(动轮，并根据式(10)，得到参与啮合的主动轮的 基体刚度kfp(τ)和参与啮合的从动轮的基体刚度kfg(τ)：
用来计算轮齿刚度的up和ug是施加之后求得的位移；是用来求齿 基刚度的ufp和ufg是施加之后求得的位移。其位移的方向都是沿着啮合 线方向，只是数值不同。
步骤43：确定齿轮副的总啮合刚度：
根据式(11)，确定齿轮副的总啮合刚度k(τ)
步骤5：确定齿根弯曲应力
参照步骤42中的方法约束内孔节点，根据载荷分配系数在啮合点施 加啮合力，进行静力学求解后获得各个单元的节点位移δe。
进一步地，根据式(12)，得到单个单元的应变ε：
ε＝Bδe (12)
其中δe为各单元节点的位移列阵，B为应变矩阵。
需要说明的是，δe是分别针对各个有限元单元的，为4*1的位移列阵。 而第一位移列阵、第二位移列阵是针对整体有限元模型的，为2n*1的列 阵，其中，n为该本地有限元模型中的总节点数。
根据式(13)，得到单个单元的应力σ：
其中，D为弹性矩阵；应力σ具有σx，σy，τxy这三个分量。
进一步地，采用von Mises应力准则来确定弯曲应力。
具体地，将不同单元的同一节点的应力作平均化处理，利用材料力 学基本理论求得各节点处的von Mises应力。
优选地，提取主动齿轮受拉侧过渡曲线上应力最大值作为齿根弯曲 应力，记为σb0(0上。也即主动齿轮受拉侧过渡曲线上的最大应力。
需要说明的是，针对每一个参与啮合的轮齿，其参与啮合一侧齿面 是受拉侧；其不参与啮合的一侧是受压侧。
进一步地，根据式(14)引入修正系数λ对靠近齿根部位的应力进行 修正：
式(14)中τ为无量纲啮合时间，其上限值为齿轮的重合度，也即无 量纲啮合周期T。
最后，根据式(15)，确定齿根弯曲应力为：
σb＝λσb0 (15)
以下从文献(Fernández A,Iglesias M,De-Juan A,et al.Gear transmission dynamic:Effects of tooth profile deviations and support flexibility[J].Applied Acoustics,2014,77(3):138-149。)中提取算例，验证 本发明实施例的齿轮副啮合特性分析方法。
表1中列出了这对具有较大修形量的直齿轮副的基本参数，其中， 修形量CRT＝15μm，修形长度ΔLT＝3mm。
分别利用ANSYS有限元建模及分析软件和本发明提出的方法分析 该齿轮副的时变啮合刚度和弯曲应力。
表1齿轮副齿轮参数
图6示出了在一组等量增幅的扭矩下这对具有较大修形量的直齿轮 副的时变啮合刚度。如图6所示，随着扭矩持续增大，在一个啮合周期 内，ANSYS有限元建模及分析方法和本发明实施例的方法尽管在两侧的 部分边缘点处存在不一致的情形，但在决定分析精度的区域，变化趋势 和数值量级均基本无差别。因此，本发明实施例的方法分析齿轮副啮合 刚度的精度能满足使用需求。
另一方面，本发明实施例的方法消耗的计算时间为3分钟，而ANSYS 接触有限元方法的计算时间为7.5小时。本方法具有和ANSYS接触有限 元相当的精度，但效率远高于ANSYS接触有限元方法。
当齿轮的各轮齿的均是健康的，载荷分配系数大致为中心对称。当 载荷作用点逐渐向齿顶移动时，力臂逐渐增大，弯曲应力增大，但是过 了单齿啮合区上界点之后，进入了双齿区，分配到的载荷小了，因此弯 曲应力有所减小。
如图7所示，在这对直齿轮在不修形时，本发明实施例的方法确定 的载荷分配系数和弯曲应力在趋势和数值上均与ANSYS接触有限元方 法基本保持一致。因此，本发明实施例的方法分析弯曲应力的精度能满 足使用需求。
图8中则示出了这对直齿轮在不修形时，主动轮的弯曲应力云图。 在啮合位置靠近齿根处，两种方法解算得到的弯曲应力有较大的不同； 随着啮合位置接近齿顶，两种方法解算得到的齿根弯曲应力趋于一致。 图8表明，本发明提出的向啮合点施加力载荷并求解线性方程的方式， 精度接近ANSYS接触有限元方法。
另一方面，完成一个啮合周期T内的应力分析，本方法耗时2分钟， 而ANSYS接触有限元分析方法的计算时间为2小时。本方法具有和 ANSYS接触有限元相当的精度，但效率远高于ANSYS接触有限元方法。
综上，本发明实施例的方法与ANSYS有限元方法相比，时变啮合刚 度的计算误差较小，其精度可以满足一般情况下对啮合刚度的精度要求； 弯曲应力的计算误差较小，其精度可以满足一般情况下对弯曲应力的精 度要求。而本发明实施例的方法的计算时间却远远小于ANSYS有限元方 法，因此本发明实施例的方法在满足精度要求的同时具有较高的计算效 率。
最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方 案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明， 本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技 术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这 些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方 案的范围。