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

Abstract

本发明涉及一种考虑复杂基体与裂纹扩展的直齿轮啮合刚度计算方法，该方法基于有限元理论获得啮合齿轮主从动轮的总体柔度矩阵，确定每个啮合位置下可能接触点的整体柔度矩阵；引入非线性赫兹接触理论，计算每个啮合位置下可能接触点的接触柔度矩阵；将可能接触点的整体柔度矩阵、接触柔度矩阵和初始间隙向量引入变形协调方程，计算该啮合位置的啮合刚度。该方法可以同时考虑由断裂力学获得的裂纹扩展路径，且可同时考虑复杂基体结构(包括腹板结构和减重孔结构)对直齿轮啮合刚度的影响。本发明方法可以计算同时考虑齿轮腹板结构、减重孔结构和裂纹扩展路径的直齿轮啮合刚度。本发明方法结果采用三维接触有限元方法进行了验证，结果说明本发明方法对于计算含裂纹的复杂基体齿轮啮合刚度存在较高精度。

Description

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

技术领域

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

背景技术

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

1.基于有限元分析软件

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

2.势能法

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

3.实验法

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

发明内容

(一)要解决的技术问题

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

(二)技术方案

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

U＝K^-1F (2)

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

u＝u_xcosα_n+u_ysinα_n (3)

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

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

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

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

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

-(λ_c+λ_b)F+x_s＝ε (5)

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

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

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

λ_c＝diag(λ_c1,λ_c2,λ_c3...λ_ci...λ_cn) (7)

其中，λ_ci表示在可能接触点i的接触柔度；F_i为可能接触点i的接触 力；E，L分别表示弹性模量和齿宽；

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

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

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

其中，ε＝[ε₁,ε₂,…ε_i,…ε_n]是初始间隙向量；

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

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

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

S402、根据步骤S401获得的静态传递误差x_s，计算在接触位置的直 齿轮啮合刚度，按公式(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软件进行应力强度因子K_I和K_II的求解。基 于ANSYS中获得的应力强度因子K_I和K_II，根据最大切向应力断裂准则 进行裂纹扩展角θ_n的计算，其计算方法为：

步骤2.3：根据初始裂纹参数、假定裂纹扩展增量与步骤2.2中裂纹 扩展角，进行裂纹扩展模拟，裂纹扩展模拟过程见图2c。裂纹扩展角θ_n和裂纹扩展增量l_s共同决定了裂纹扩展方向与裂纹平滑程度。根据裂纹 初始参数(初始裂纹位置角θ，初始裂纹长度l_I、初始裂纹方向α)在裂纹 初始位置建立初始裂纹AB(见图2b)，并在裂纹尖端B点建立局部平面坐标系(该坐标系x轴沿裂纹AB方向，y轴沿裂纹AB垂直方向)。根据裂 纹扩展角θ_n与裂纹扩展增量l_s，可以确定下一步裂纹扩展步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软件中，对于主从动轮轴孔上节点的自由度进行全约束，并对 主从动轮上沿啮合方向的接触点施加单位力F_u，见图3(a)。

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

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

U＝K^-1F (2)

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

u＝u_xcosα_n+u_ysinα_n (3)

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

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

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

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

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

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

-(λ_c+λ_b)F+x_s＝ε (5)

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

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

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

λ_c＝diag(λ_c1,λ_c2,λ_c3...λ_ci...λ_cn) (7)

其中，λ_ci表示在可能接触点i的接触柔度；F_i为可能接触点i的接触 力；E，L分别表示弹性模量和齿宽。

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

其中，F是总的接触力(F＝T/r_b1，T是扭矩，r_b1是主动轮基圆半径)； F_i表示可能接触点i的接触力。

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

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

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

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

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

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

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

实施例2

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

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

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

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

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

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

步骤5.1：采用表1齿轮参数，用本发明提出方法，现有技术中的能 量法和有限元方法(参见Ma H,Zeng J,Feng R.,et al.An improved analytical method for meshstiffness calculation of spur gears with tip relief [J].Mechanism and MachineTheory,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、减重孔位置半径r_sp和减重孔个数分别为20mm、 60mm和5。本发明方法算获得的刚度结果变化趋势与有限元结果相似，其相对有限元方法的的最大误差为3.5％。同时可知，减重孔会引起齿轮 刚度存在一个整体的波动。

表3齿轮副基本参数

表4不同案例的腹板参数

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

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

本发明扩展获得的裂纹路径采用扩展有限元方法和实验方法进行了 验证(试验台采用文献Li M,Xie L,Ding L.Load sharing analysis and reliability predictionfor planetary gear train of helicopter[J].Mechanism and Machine Theory,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复杂基体齿轮副参数

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

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

Claims (6)

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

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

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

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

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

2.根据权利要求1所述的直齿轮啮合刚度计算方法，其特征在于，在步骤S1中，所述裂纹直齿轮副的基本参数包括主从动轮齿数、模数、齿宽、分度圆压力角、弹性模量、泊松比、内孔径；所述腹板参数包括腹板內圆半径、腹板外圆半径、腹板厚度；所述初始裂纹参数包括初始裂纹位置角θ、初始裂纹长度l_I和初始裂纹方向α；减重孔参数包括圆形减重孔参数和形减重孔参数，其中，圆形减重孔参数包括：减重孔直径Φ_S、减重孔位置半径r_sp和圆形减重孔个数；扇形减重孔参数包括：减重孔内径r_si、外径r_so、扇形减重孔展开角和扇形减重孔个数。

3.根据权利要求1所述的直齿轮啮合刚度计算方法，其特征在于，所述步骤S2具体包括如下步骤：

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

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

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

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

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

4.根据权利要求1所述的直齿轮啮合刚度计算方法，其特征在于，在步骤S203中，所述裂纹扩展模拟通过初始裂纹参数确定初始裂纹位置与扩展方向，并通过每步扩展计算获得的裂纹扩展角θ_n和裂纹扩展增量l_s来决定下一步裂纹扩展，重复该过程，完成整个裂纹扩展模拟。

5.根据权利要求1所述的直齿轮啮合刚度计算方法，其特征在于，所述步骤S3具体包括如下步骤：

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

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

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

U＝K^-1F (2)

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

u＝u_xcosα_n+u_ysinα_n (3)

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

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

其中，λ^Pr代表主动轮粗略的柔度矩阵；N表示主动轮齿面1-3上的节点总数；λ_ij表示施力节点为j时，节点i的柔度。

6.根据权利要求1所述的直齿轮啮合刚度计算方法，其特征在于，在步骤S4中，包括如下步骤：

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

-(λ_c+λ_b)F+x_s＝ε (5)

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

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

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

λ_c＝diag(λ_c1,λ_c2,λ_c3...λ_ci...λ_cn) (7)

其中，λ_ci表示在可能接触点i的接触柔度；F_i为可能接触点i的接触力；E，L分别表示弹性模量和齿宽；

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

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

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

其中，ε＝[ε₁,ε₂,…ε_i,…ε_n]是初始间隙向量；

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

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

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

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