考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法与流程

本发明涉及一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，属于机械动力学技术领域。

背景技术：

目前，现有的可以同时考虑复杂齿基与裂纹扩展路径的直齿轮啮合刚度计算方法主要有以下几种方法：

1.基于有限元分析软件

在有限元软件ansys，abaqus等中建立含复杂基体和裂纹扩展路径的三维有限元模型，并在软件中，选择合适的单元及材料参数，完成三维模型网格划分，建立合适的接触对，设置合适的约束并选择适当的求解方法对含复杂基体与裂纹路径的直齿轮啮合刚度进行计算。利用现有的有限元分析软件进行啮合刚度的计算，往往建模过程复杂且繁重，计算效率很低。

2.势能法

基于前人研究的不断发展与完善，势能法逐渐趋于完善。对于势能法，学者们将齿轮对的时变啮合刚度主要分为接触刚度、弯曲刚度、剪切刚度、轴向压缩刚度，以及齿轮基体刚度几大部分。通过各分支刚度之间的串并联关系，计算直齿轮副的啮合刚度。势能法计算啮合刚度具有比较高的计算效率，但是相对有限元方法，其在基体刚度与接触刚度的计算中存在一定误差。同时，其很难准确考虑复杂齿轮基体和裂纹扩展路径对啮合刚度的影响。

3.实验法

一些学者采用实验方法进行齿轮啮合刚度的计算。实验方法主流的方法有：光弹测量技术、应变仪技术、激光位移传感器和数字图像关联技术等。实验方法可以更加准确的考虑齿轮不对中、齿面误差等真实因素对啮合刚度的影响。但是实验方法存在操作难度大，经济投入高，干扰因素多等问题，因此，实验测量的结果往往存在一定的误差。

技术实现要素：

(一)要解决的技术问题

为了解决现有技术的上述问题，本发明提供一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，该方法在保证计算精度的前提下，充分考虑齿轮复杂基体结构(腹板和减重孔)、裂纹扩展路径(轮齿裂纹和基体裂纹)等对直齿轮啮合刚度的影响；其在保证较高计算效率的同时保证了较高的计算精度。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，包括如下步骤：

s1、确定直齿轮副参数，参数包括基本参数、腹板参数、减重孔参数和初始裂纹参数；

s2、根据所述直齿轮副的基本参数与初始裂纹参数，采用有限元软件进行裂纹扩展路径的模拟；

s3、基于直齿轮副参数，借助有限元软件建立考虑减重孔参数的有限元网格模型，将含裂纹路径或无裂纹路径和减重孔的网格模型节点与单元信息导入matlab软件，基于有限元理论即通过控制有限元理论中平面四节点单元厚度考虑齿轮辐板参数与承载轮齿接触分析方法，进行包含复杂基体与裂纹路径的主从动齿轮柔度矩阵的计算；

s4、根据s3中获得的主从动轮柔度矩阵，将其引入变形协调方程，并结合赫兹接触理论，进行静态传递误差xs的求解，并进行啮合刚度k的计算。

如上所述的直齿轮啮合刚度计算方法，优选地，在步骤s1中，所述裂纹直齿轮副的参数包括基本参数：主从动轮齿数、模数、齿宽、分度圆压力角、弹性模量、泊松比、内孔径；腹板参数：腹板內圆半径、腹板外圆半径、腹板厚度；初始裂纹参数：初始裂纹位置角θ、初始裂纹长度li和初始裂纹方向α；减重孔参数包括圆形减重孔参数和形减重孔参数，其中，圆形减重孔参数包括：减重孔直径φs、减重孔位置半径rsp和圆形减重孔个数；扇形减重孔参数包括：减重孔内径rsi、外径rso、扇形减重孔展开角和扇形减重孔个数。

如上所述的直齿轮啮合刚度计算方法，优选地，所述步骤s2具体包括如下步骤:

s201、根据直齿轮副的基本参数与初始裂纹参数，在ansys软件中建立含初始裂纹的有限元网格模型；

s202、采用ansys软件进行应力强度因子ki和kii的求解；基于ansys中获得的应力强度因子ki和kii，根据最大切向应力断裂准则进行裂纹扩展角θn的计算，其计算方法为：

s203、根据初始裂纹参数、假定裂纹扩展增量与步骤s202中裂纹扩展角，进行裂纹扩展模拟；

s204、当步骤s203中完成单步的裂纹扩展时，建立新裂纹长度下的含奇异单元的有限元模型，依次重复步骤s202-s203，完成下一步裂纹扩展过程的模拟；

s205、重复上述裂纹扩展过程，完成裂纹扩展路径的模拟。

如上所述的直齿轮啮合刚度计算方法，优选地，在步骤s203中，所述裂纹扩展模拟主要通过初始裂纹参数确定初始裂纹位置与扩展方向，并通过每步扩展计算获得的裂纹扩展角θn和裂纹扩展增量ls来决定下一步裂纹扩展，重复该过程，完成整个裂纹扩展模拟。

如上所述的直齿轮啮合刚度计算方法，优选地，所述步骤s3具体包括如下步骤：

s301、利用s2中获得的裂纹扩展路径和s1中的减重孔参数，在ansys软件中建立包含裂纹路径和减重孔结构的有限元网格模型，将该网格模型的节点与单元信息导入matlab软件；

s302、在所述matlab软件中，通过有限元理论控制平面单元的厚度，实现辐板结构的考虑，进而，利用有限元理论，对主从动轮轴孔上节点的自由度进行全约束，并对主从动轮上沿啮合方向的接触点施加单位力fu；

s303、根据边界条件与有限元理论，按公式(2)求解全局变形矢量u，

u＝k^-1f(2)

其中，f是整体节点力矢量；k是局部刚性区与内孔约束的整体刚度矩阵，根据全局变形矢量u，主动轮轮齿各齿面上每个节点的总变形可以获得：

u＝uxcosαn+uysinαn(3)

其中，u表示提取变形节点的总变形；ux,uy分别表示力施加接触节点时，提取位移节点沿水平和竖直方向的变形；αn是啮合压力角；

当主动轮轮齿齿面上所有节点依次施加一次单位力，获得主动轮粗略的柔度矩阵λ^pr，其表达如公式(4)所示：

其中，λ^pr代表主动轮粗略的柔度矩阵；n表示主动轮齿面上的节点总数；λij表示施力节点为j时，节点i的柔度。

如上所述的直齿轮啮合刚度计算方法，优选地，在步骤s4中，包括如下步骤：

s401、进行静态传递误差xs的求解，先引入变形协调方程，如公式(5)所示：

-(λc+λb)f+xs＝ε(5)

其中，xs表示静态传递误差；ε表示初始间隙向量，λb表示所有可能接触点的总体柔度矩阵，其如公式(6)所示：

其中，n表示可能接触点总数；表示施力点为j时，主动轮提取点i的柔度；表示在该相应接触位置的从动轮节点柔度，其中和都从总体柔度矩阵λ^p和λ^g中获取；

当轮齿处于n对轮齿接触时，接触柔度矩阵如公式(7)所示为：

λc＝diag(λc1,λc2,λc3...λci...λcn)(7)

其中，λci表示在可能接触点i的接触柔度；fi为可能接触点i的接触力；e，l分别表示弹性模量和齿宽；

f表示n齿对接触时，所有可能接触点构成的法向接触力矢量，其如公式(9)所示：

其中，f是总的接触力；fi表示可能接触点i的接触力；

对于健康齿轮，基于公式(5)-(9)，变形协调方程转换为如公式(10)所示：

其中，ε＝[ε1,ε2,…εi,…εn]是初始间隙向量；

为了求解公式(10)，将给定接触力初值的方式进行迭代，接触力的初值设定为：

在迭代过程中，如果fi<0，则接触点i被认为是虚假接触，则公式(10)变形为如公式(12)所示：

通过上式继续进行迭代，直至f1,f2,…,fn>0且f1,f2,…,fn保持不变，完成迭代过程；最终，通过变形协调方程的求解可以获得所有可能接触点的接触力矢量f与静态传递误差xs；

s402、根据步骤s401获得的静态传递误差xs，计算在接触位置的直齿轮啮合刚度，按公式(13)计算；

需要进行说明，本发明方法可以考虑腹板结构(需要用到腹板参数)，也可以考虑减重孔结构(减重孔参数)，当然这些也不是必要的齿轮参数。也就说，本方法可以不考虑腹板结构、减重孔结构和裂纹，同样可以计算健康实心齿轮副的啮合刚度。因此，在后边具体实施方法中，研究这几个因素(腹板结构、减重孔结构和裂纹)影响的时候，图7是不考虑腹板结构和减重孔结构的健康实心齿轮啮合刚度结构；图9a和9b是考虑腹板结构的健康齿轮刚度结果；图10是考虑减重孔结构的健康齿轮啮合刚度结构；图12为不考虑腹板结构和减重孔结构，含裂纹的直齿轮啮合刚度结构。图15,16为含裂纹，腹板结构，减重孔结构计算获得的刚度结构。

(三)有益效果

本发明的有益效果是：

本发明提供一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，其是基于有限元理论与受载轮齿接触分析方法建立了直齿轮啮合刚度的计算方法，考虑了采用断裂力学方法扩展获得的较为真实的裂纹扩展路径，非假设的裂纹路径，可计算考虑真实裂纹扩展路径的直齿轮啮合刚度。对于单位集中力载荷在齿面造成的不合理的局部变形，本发明采用在啮合位置建立半径为0.2m(m为齿轮模数)，弹性模量调大1000倍的局部刚性区来防止啮合位置施加集中力造成的不合理的局部变形。通过控制有限元理论中平面单元的厚度来实现齿轮腹板结构的模拟。该方法对于模拟齿轮腹板结构具有一定精度，该结果已通过三维接触有限元方法进行了验证，通过控制有限元理论中平面单元的厚度与基体结构中单元划分，本发明方法可以模拟齿轮腹板结构与减重孔共同构成的复杂齿基结构，三维接触有限元方法进行验证，结果说明本发明方法对于计算含裂纹的复杂基体齿轮啮合刚度存在较高精度。

本发明相对于背景技术中提到的有限元软件计算方法，具有更高的计算效率，且精度较高；本发明提出方法相对于势能法具有更高的精度，其结果已经在图7证明。同时，本发明方法相对于实验法，经济成本耗费低，且效率高。综上，本发明的方法具有较高的计算精度与效率。

附图说明

图1为本发明具体实施方式中的考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算的流程图；

图2为含奇异单元的裂纹扩展有限元模型与裂纹扩展过程示意图，其中，(a)裂纹齿轮有限元模型示意图，(b)裂纹扩展过程示意图；(c)裂纹扩展模拟过程；

图3为本文方法模型边界条件示意图，其中，(a)整体约束示意图，(b)齿面载荷示意图；

图4为本文方法与有限元模型复杂齿基示意图，其中，(a)本文复杂齿基示意图，(b)有限元模型示意图；

图5为本文方法主动轮轮齿1-3齿面柔度矩阵示意图；

图6为考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法流程图；

图7为本发明方法、势能法和有限元方法的啮合刚度对比；

图8为腹板齿轮结构与减重孔齿轮结构示意图，其中，(a)腹板齿轮基体结构参数示意图，(b)腹板齿轮三维有限元模型示意图，(c)减重孔齿轮基体结构参数示意图，(d)减重孔齿轮三维有限元模型示意图；

图9为齿轮不同腹板参数下的啮合刚度，其中，(a)本发明方法的刚度结果，(b)有限元方法结果；

图10为含圆形减重孔的齿轮啮合刚度；

图11为实心基体齿轮不同裂纹长度下裂纹路径示意图，其中，(a)本发明方法中获得裂纹路径，(b)扩展有限元方法结果，(c)实验方法结果；

图12为实心基体齿轮不同裂纹长度下刚度结果，其中，(a)本发明方法的刚度结果，(b)有限元方法结果；

图13为扇形齿轮基体结构的示意图，其中，(a)扇形齿轮基体结构参数示意图，(b)扇形基体齿轮三维有限元模型示意图；

图14为不同齿基位置进行裂纹扩展获得的裂纹扩展路径，其中，(a)基体裂纹路径，(b)轮齿裂纹路径；

图15为基体裂纹齿轮啮合刚度结果，其中，(a)本发明方法获得的基体裂纹齿轮啮合刚度结果，(b)有限元方法获得的基体裂纹齿轮啮合刚度结果；

图16为轮齿裂纹齿轮啮合刚度结果，(a)本发明方法获得的轮齿裂纹齿轮啮合刚度结果，(b)有限元方法获得的轮齿裂纹啮合刚度结果。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

实施例1

一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，如图1所示，主要包括如下步骤：步骤1：确定直齿轮副的基本参数、腹板参数、减重孔参数和初始裂纹参数；

步骤2：确定初始裂纹参数，进行裂纹扩展路径的模拟。具体采用如下步骤：

步骤2.1：根据步骤1中确定的直齿轮副参数，在ansys软件中建立含初始裂纹的有限元模型，建立的包含奇异单元的有限元模型如图2a所示；

步骤2.2：采用ansys软件进行应力强度因子ki和kii的求解。基于ansys中获得的应力强度因子ki和kii，根据最大切向应力断裂准则进行裂纹扩展角θn的计算，其计算方法为：

步骤2.3：根据初始裂纹参数、假定裂纹扩展增量与步骤2.2中裂纹扩展角，进行裂纹扩展模拟，裂纹扩展模拟过程见图2c。裂纹扩展角θn和裂纹扩展增量ls共同决定了裂纹扩展方向与裂纹平滑程度。根据裂纹初始参数(初始裂纹位置角θ，初始裂纹长度li、初始裂纹方向α)在裂纹初始位置建立初始裂纹ab(见图2b)，并在裂纹尖端b点建立局部平面坐标系(该坐标系x轴沿裂纹ab方向，y轴沿裂纹ab垂直方向)。根据裂纹扩展角θn与裂纹扩展增量ls，可以确定下一步裂纹扩展步bc(其中裂纹扩展角θn是x轴与下一步裂纹扩展方向之间的夹角，见图2(b))。

步骤2.4：当步骤2.3中完成单步的裂纹扩展时，建立新裂纹长度下的含奇异单元的有限元模型。依次重复步骤2.2-2.3，完成下一步裂纹扩展过程的模拟。

步骤2.5：重复上述裂纹扩展过程，完成裂纹扩展路径的模拟。

步骤3：基于直齿轮副参数，借助有限元软件建立考虑减重孔的平面网格模型，将含裂纹路径(裂纹路径由步骤2获得)或无裂纹路径和减重孔的网格模型节点与单元信息导入matlab软件，基于有限元理论与承载轮齿接触分析方法，进行包含复杂基体的主从动齿轮柔度矩阵的计算。

步骤3.1：根据ansys软件中建立的包含裂纹扩展路径(或无裂纹路径)和减重孔的有限元网格模型，将该网格模型的节点与单元信息导入matlab软件。

步骤3.2：基于有限元理论，进行考虑复杂齿基的总刚矩阵组集。在matlab软件中，对于主从动轮轴孔上节点的自由度进行全约束，并对主从动轮上沿啮合方向的接触点施加单位力fu，见图3(a)。

需要注意的是，如果在轮齿上施加单位集中力，此时会在齿面接触点形成较大的局部变形，该变形不能真正表征齿轮的整体变形。因此，本发明提出了在施加力位置建立刚度区的方式来解决集中力导致的局部变形。根据大量数据的验证，确定取刚性区半径为0.2m(m为齿轮模数)的圆形区域，将该位置内单元的弹性模量扩大1000倍，即，认为是刚性的，如图3(b)所示。以主动轮的轮齿2为例，在其齿面有多个啮合节点。假设在如图3(b)所示位置啮合，在主从动轮的该位置分别建立半径为0.2m的圆形刚性区，并在主从动轮该啮合位置施加一个沿啮合线方向的单位力fu。根据有限元理论与边界条件，进行平面四节点单元刚度矩阵组集，获得考虑局部刚性区与内孔约束的整体刚度矩阵k。需要说明的是，本发明对于二维齿轮采用了带厚度的平面四节点单元进行齿轮结构的模拟。对于健康实心齿轮，采用单元的厚度模拟了齿宽。为了模拟齿轮的腹板厚度，同样采用修改单元厚度的方式来模拟。如图4(a)所示，通过解析中节点坐标选中腹板位置的单元(见图4(a)中深色单元，进行腹板位置的单元厚度控制，从而实现真实直齿轮腹板结构的模拟。其真实的直齿轮腹板结构可以参考图4(b)。

步骤3.3：根据边界条件与有限元理论，求解全局变形矢量u:

u＝k^-1f(2)

其中，f是整体节点力矢量；k是考虑局部刚性区与内孔约束的整体刚度矩阵(根据有限元理论单元刚度矩阵组集获得)；为了提高计算效率，编程时采用了cholesky分解与稀疏矩阵方法求解。根据全局变形矢量u，主动轮轮齿1-3齿面上每个节点的总变形可以获得：

u＝uxcosαn+uysinαn(3)

其中，u表示提取变形节点的总变形；ux,uy分别表示力施加接触节点时，提取位移节点沿水平和竖直方向的变形；αn是啮合压力角。

当单位力施加在轮齿2的接触节点，在主动轮轮齿1-3中齿面每个节点的总变形都可以获得(参考图5)。由于在齿面施加的是单位力，因此，节点的变形也可以被称为节点的柔度。进而，当主动轮轮齿1-3齿面上所有节点依次施加一次单位力，可以获得主动轮粗略的柔度矩阵λ^pr，其可以表达为：

其中，λ^pr代表主动轮粗略的柔度矩阵；n表示主动轮齿面1-3上的节点总数；λij表示施力节点为j时，节点i的柔度，如图5所示。

由于齿轮齿面的节点个数有限，因此很难准确获得每个啮合位置的柔度结果。因此，本发明采用了插值的方法，获得了更为准确的齿面柔度结果(见图5)。插值之后，主动轮的柔度矩阵为λ^p。采用相似的方法可以获得从动轮的柔度矩阵λ^g。

步骤4：根据变形协调方程和赫兹接触理论，进行静态传递误差xs的求解，并进行啮合刚度k的计算。

步骤4.1：进行静态传递误差xs的求解。首先引入变形协调方程，其表达为：

-(λc+λb)f+xs＝ε(5)

其中，xs表示静态传递误差；ε表示初始间隙向量，其与装配误差、加工误差、修形量及理论分离距离等相关；λb表示所有可能接触点的总体柔度矩阵，其表示为：

其中，n表示可能接触点总数；表示施力点为j时，主动轮提取点i的柔度；表示在该相应接触位置的从动轮节点柔度，其中和都可以从总体柔度矩阵λ^p和λ^g中获取。

当轮齿处于n对轮齿接触时，接触柔度矩阵可以表示为：

λc＝diag(λc1,λc2,λc3...λci...λcn)(7)

其中，λci表示在可能接触点i的接触柔度；fi为可能接触点i的接触力；e，l分别表示弹性模量和齿宽。

f表示n齿对接触时，所有可能接触点构成的法向接触力矢量，其可以表示为：

其中，f是总的接触力(f＝t/rb1，t是扭矩，rb1是主动轮基圆半径)；fi表示可能接触点i的接触力。

对于健康齿轮，基于式(5)-(9)，变形协调方程可以转换为：

其中，ε＝[ε1,ε2,…εi,…εn]是初始间隙向量。一般而言，直齿轮主要存在单齿接触与双齿接触。因此，可能的接触点数n为1或2。然而，本发明中为了考虑延长啮合的影响，假设可能接触点数n为3，并在初始接触间隙ε中考虑了延迟啮合与提前啮入的分离距离sa和sr。

为了求解式(10)，本发明将给定接触力初值的方式进行迭代，接触力的初值设定为：

在迭代过程中，如果fi<0，则接触点i被认为是虚假接触(接触点i不处于接触状态)。则式(10)可以变形为：

通过上式继续进行迭代，直至f1,f2,…,fn>0且f1,f2,…,fn保持不变，完成迭代过程。最终，通过变形协调方程的求解可以获得所有可能接触点的接触力矢量f与静态传递误差xs。

步骤4.2：基于获得的静态传递误差xs，进行直齿轮啮合刚度的计算。因此，在接触位置的啮合刚度可以表示为：

综上，考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法的流程图如图6所示。

实施例2

一种考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，本实施例是在实施例1的基础上，包括以下步骤：

步骤1：确定直齿轮副的基本参数及腹板位置和减重孔参数；

步骤2：确定初始裂纹参数，进行裂纹扩展路径的模拟。具体参考发明内容步骤2；

步骤3：基于有限元理论与承载轮齿接触分析方法，进行包含复杂基体与裂纹的主从动齿轮柔度矩阵的计算。具体实施方法参考发明内容中步骤3；

步骤4：采用变形协调方程和赫兹接触理论，进行静态传递误差xs的求解，进而计算完成考虑复杂基体与裂纹扩展路径的齿轮啮合刚度k。具体实施方法参考发明内容中步骤4；

步骤5：基于步骤1-4中介绍的考虑复杂基体与裂纹扩展路径的直齿轮啮合刚度计算方法，本步骤主要采用有限元方法进行本发明方法的验证，并对复杂基体与裂纹参数进行影响分析。

步骤5.1：采用表1齿轮参数，用本发明提出方法，现有技术中的能量法和有限元方法(参见mah,zengj,fengr.,etal.animprovedanalyticalmethodformeshstiffnesscalculationofspurgearswithtiprelief[j].mechanismandmachinetheory,2016,98:64-80.6)的获得的刚度结果，如图7所示。采用能量法和有限元方法来验证本发明提出方法对计算健康直齿轮啮合刚度的准确性。通过比较可知，三种方法刚度变化趋势相似。同时，双齿啮合区和单齿啮合区中间时刻a，b时刻的刚度及其误差见表2。由表格数据知，在a、b时刻，提出方法相对于有限元方法结果的最大误差为1.023％，其误差值相对能量法较小。

表1齿轮副基本参数

表2在a、b时刻啮合刚度的误差

注：表中的数据均为其他方法相对于有限元方法结果的误差。

步骤5.2：采用本发明方法与有限元方法分析了腹板结构参数(腹板厚度和轮缘厚度)和圆形减重孔对齿轮啮合刚度的影响。腹板结构齿轮基体结构与减重孔齿轮基体结构示意图如图8所示，采用的齿轮参数见表3，腹板参数见表4，齿轮不同腹板参数下的啮合刚度如图9所示；通过案例2和例3的比较，轮缘厚度降低10mm，导致双齿区啮合刚度降低3.32％；而通过案例3和案例4的比较，腹板厚度降低2mm，导致双齿区啮合刚度下降7.67％。因此，可知腹板厚度对啮合刚度的影响比轮缘厚度的影响更明显。含圆形减重孔的齿轮啮合刚度如图10所示，其采用的圆形减重孔直径φs、减重孔位置半径rsp和减重孔个数分别为20mm、60mm和5。本发明方法算获得的刚度结果变化趋势与有限元结果相似，其相对有限元方法的的最大误差为3.5％。同时可知，减重孔会引起齿轮刚度存在一个整体的波动。

表3齿轮副基本参数

表4不同案例的腹板参数

步骤5.3：基于步骤2中获得的裂纹扩展路径，采用本发明方法与有限元方法分析了裂纹扩展路径对实心基体齿轮啮合刚度的影响。

具体采用表1的齿轮参数，以及初始裂纹参数(为了保证仿真方法的裂纹初始参数与实验的初始裂纹参数一致，因此，本发明方法与扩展有限元方法模拟该条件下的裂纹扩展路径时，裂纹初始位置选择为健康齿轮齿根应力最大位置，且其初始裂纹长度li与初始裂纹方向角α分别为0.2mm和45°)，实心基体齿轮不同裂纹长度下裂纹路径如图11所示。

本发明扩展获得的裂纹路径采用扩展有限元方法和实验方法进行了验证(试验台采用文献lim,xiel,dingl.loadsharinganalysisandreliabilitypredictionforplanetarygeartrainofhelicopter[j].mechanismandmachinetheory,2017,115:97-113中的方法)。采用本发明提出的啮合刚度计算方法，计算了不同裂纹长度4.0mm、8.0mm和12.0mm的齿轮啮合刚度。本发明方法与有限元方法在不同裂纹长度下的刚度结果如图12所示，其中，(a)为本发明方法的刚度结果，(b)为有限元方法结果。两种方法获得的刚度计算结果变化趋势相似。且随着裂纹长度的增加，齿轮啮合刚度逐渐减小，同时可以看出主动轮裂纹对双齿区ii的影响比双齿区i刚度影响较大。这是由于当处于主动轮双齿区ii的时候，啮合位置接近裂纹齿的齿顶，啮合力的作用力臂较长，会造成裂纹齿较大的变形，进而造成啮合刚度减小。

步骤5.4：本部分主要对于不同裂纹路径的复杂齿基直齿轮啮合刚度进行计算分析。采用本发明方法与有限元方法分析了复杂基体结构对裂纹扩展路径的影响，以及裂纹扩展路径对复杂齿体齿轮啮合刚度的影响。采用表5的齿轮、辐板和减重孔参数，建立如图13所示从动轮扇形减重孔和腹板结合的基体结构形式，基于该齿轮基体结构，进行裂纹路径扩展。其中，图13(a)为扇形齿轮基体结构参数示意图，图13(b)为扇形基体齿轮三维有限元模型示意图。采用裂纹初始参数(见表5)：初始裂纹位置角、初始裂纹长度和初始裂纹扩展角分别为：104°、0.3mm、45°，在不同的基体位置进行了裂纹扩展，获得不同的裂纹路径如图14所示，(a)基体裂纹路径，(b)轮齿裂纹路径。对于裂纹类型i，如图14(a)，该裂纹类型裂纹长度分别为0,4.0,8.0mm(见图14a)下的齿轮啮合刚度如图15所示，其中(a)为本发明方法获得的基体裂纹齿轮啮合刚度结果，(b)为有限元方法获得的基体裂纹齿轮啮合刚度结果。对于裂纹类型ii，由于腹板结构的支撑作用，裂纹主要为轮齿裂纹。对于该裂纹类型裂纹长度分别为0、4.0、8.0、12.0mm的啮合刚度进行了计算，结果如图16所示，其中(a)为本发明方法获得的轮齿裂纹齿轮啮合刚度结果，(b)为有限元方法获得的轮齿裂纹啮合刚度结果。

对于本发明方法计算的结果，采用了有限元方法进行了验证。通过本发明方法与有限元方法的结果比较，两种方法的啮合刚度变化趋势相同，且本发明方法的最大计算误差仅为4％，这证明了本发明方法对于计算含裂纹的复杂基体齿轮啮合刚度存在较高精度。

表5复杂基体齿轮副参数

说明：齿顶高系数ha^*，顶隙系数c^*分别为1.00和0.25。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明做其它形式的限制，任何本领域技术人员可以利用上述公开的技术内容加以变更或改型为等同变化的等效实施例。但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。