基于解析-有限元的复杂齿基直齿圆柱齿轮啮合刚度计算方法与流程

本发明属于齿轮精密传动技术和力学分析领域，涉及基于解析-有限元的复杂齿基直齿圆柱齿轮啮合刚度计算方法，可用于直齿圆柱齿轮的动态性能分析及优化设计。

背景技术：

在齿轮传统系统中，主要包括直齿轮传动、斜齿轮传动和蜗轮蜗杆传动等多种传动方式，其中的直齿轮传动又可进一步分为直齿圆柱齿轮传动、直齿圆锥齿轮传动和直齿端面齿轮传动等形式，而直齿圆柱齿轮传动是在传动领域中运用最为广泛的传动形式之一，这种传动结构不仅具有传动平稳，冲击、振动和噪声较小等特点，还具备体积小、重量轻、传递转矩大、起动平稳和传动比分级精细的优点，因此广泛应用于高端数控装备、自动化精密机械和航空航天等领域。

近几年来，高端数控装备正不断向着高速、高精度和长寿命的方向发展，并对直齿圆柱齿轮这一关键部件的动力学特性提出了更高的要求。而啮合刚度周期性的变化是传动系统的主要激励形式之一，这种啮合刚度周期性的变化称之为时变啮合刚度，其直接影响着传动过程中的动力学特性。时变啮合刚度主要是由于直齿圆柱齿轮的重合度一般不为整数，即直齿圆柱齿轮同时参与啮合的齿数随时间呈周期性变化造成的。在直齿圆柱齿轮的实际加工中，为了获得良好的机械性能，直齿圆柱齿轮的齿基往往被加工成腹板带孔的结构。因此能否精确求解这种含有复杂齿基的直齿圆柱齿轮时变啮合刚度是研究直齿圆柱齿轮传动系统动力学特性，进行动态性能分析及优化设计的前提，所以研究含有复杂齿基的直齿圆柱齿轮时变啮合刚度的计算方法很有必要的。

从目前公开的资料来看，在计算含有复杂齿基的直齿圆柱齿轮的时变啮合刚度方面主要采用有限元法进行模拟计算。例如，2010年abbes在其发表的名称为“dynamicbehaviourmodellingofaflexiblegearsystembytheelasticfoundationtheoryinpresenceofdefects”的期刊论文中公开了一种考虑复杂齿基的直齿圆柱齿轮的动力学分析计算方法。该方法主要是采用有限元仿真的方法获得含有复杂齿基的直齿圆柱齿轮的动力学响应。但是这种方法需要建立完整的有限元模型来模拟真实的情况，需要耗费大量的计算资源，因此存在着计算效率低的不足。

技术实现要素：

本发明的目的在于克服上述现有技术存在的缺陷，提出了基于解析-有限元的复杂齿基直齿圆柱齿轮啮合刚度计算方法，旨在提高含有复杂齿基的直齿圆柱齿轮啮合刚度的计算效率和精度。

为了实现上述目的，本发明采用的技术方案包括如下步骤：

(1)基于解析法，分别计算直齿圆柱齿轮轮齿刚度，包括轮齿弯曲刚度值kb、剪切刚度值ks和径向压缩刚度值ka；

(2)基于赫兹接触理论，计算轮齿的接触刚度值kh；

(3)基于刚度串并联理论，计算单齿对啮合刚度值kt；

(4)根据含复杂齿基的直齿圆柱齿轮几何参数建立含复杂齿基直齿圆柱齿轮的几何模型，并基于有限单元法，采用solid185单元进行网格映射划分，形成含复杂齿基的直齿圆柱齿轮的有限元模型；

(5)根据齿轮啮合理论，建立刚性耦合区域；

(5.1)基于直齿圆柱齿轮的啮合理论，分别在主动齿轮和从动齿轮的啮合点位置处建立关键点，并采用mass21单元划分关键点，形成主质量点；

(5.2)基于主质量点与复杂齿基上轮齿的对应关系，采用mpc184单元将主质量点与对应的齿基部分进行刚性耦合，形成刚性耦合区域；

(6)根据齿轮啮合理论，判断主质量点对应的轮齿是否处于啮合状态，若处于啮合状态，则将处于啮合状态的轮齿对应的主质量点用combin14单元连接，并将combin14单元的刚度值设置为步骤(3)中得到的kt值；

(7)基于静态分析，计算含复杂齿基的直齿圆柱齿轮的啮合刚度值k；

(8)判断是否完成一个啮合周期的刚度计算，若没有完成一个啮合周期的刚度计算，则直齿圆柱齿轮旋转一定角度，重复步骤(1)至步骤(7)，若完成了一个啮合周期的刚度计算，则停止计算，得到含复杂齿基直齿圆柱齿轮的时变啮合刚度。

作为本发明的一种优选，步骤(1)中所述的分别计算直齿圆柱齿轮轮齿刚度，包括轮齿弯曲刚度值kb、剪切刚度值ks和径向压缩刚度值ka，其计算表达式分别为：

其中，kb为轮齿的弯曲刚度，ks为轮齿的剪切刚度，ka为轮齿的径向压缩刚度，e为弹性模量，g为剪切模量，为分度圆压力角，rb为基圆半径，τc＝αc-(θb-invαc)，为接触点处的压力角，θb为基齿角的一半，为齿顶高系数，m为模数，γ和τ为积分范围，y1为过渡曲线上任一点处的水平坐标，y2为渐开线上任一点处的水平坐标，iy1为过渡曲线上任意位置处的截面惯性矩，iy2为渐开线上任意位置处的截面惯性矩，ay1为过渡曲线上任意位置处的截面面积，ay2为渐开线上任意位置处的截面面积。

作为本发明的一种优选，步骤(2)中所述的计算轮齿的接触刚度值kh，其计算表达式为：

其中，kh为轮齿之间的接触刚度，e为弹性模量，b为切片的宽度，v为切片的泊松比。

作为本发明的一种优选，步骤(3)中所述的计算单齿对啮合刚度kt，其计算表达式为：

其中，kt为单齿对的啮合刚度，kb1为主动轮轮齿的弯曲刚度，ks1为主动轮轮齿的剪切刚度，ka1为主动轮轮齿的径向压缩刚度，kb2为从动轮轮齿的弯曲刚度，ks2为从动轮轮齿的剪切刚度，ka2为从动轮轮齿的径向压缩刚度，kh为轮齿之间的接触刚度。

作为本发明的一种优选，步骤(4)中所述的根据含复杂齿基的直齿圆柱齿轮几何参数建立含复杂齿基直齿圆柱齿轮的几何模型，其具体建立过程为：

首先，将直齿圆柱齿轮简化为一个圆柱体，圆柱体的高度等于直齿圆柱齿轮的齿宽，圆柱体的半径等于直齿圆柱齿轮的齿根圆半径；其次，根据直齿圆柱齿轮的参数，在圆柱体上完成腹板部分的建模；最后，根据直齿圆柱齿轮的参数，在腹板上完成减重孔的建模，从而形成含复杂齿基的直齿圆柱齿轮的几何模型。

作为本发明的一种优选，步骤(7)中所述的计算含复杂齿基的直齿圆柱齿轮的啮合刚度值k，其计算表达式为：

其中，t为施加在主动轮中心孔处的扭矩，rb1为主动轮中心孔半径，δθ为主动轮中心孔的变形。

本发明与现有技术相比，具有如下优点：

本发明采用解析法与有限元法相结合的方式来计算含有复杂齿基的直齿圆柱齿轮的时变啮合刚度，并且在主质量点的连接方式上采用了combin14单元进行连接，相比现有的只通过有限元法来计算含有复杂齿基的直齿圆柱齿轮的时变啮合刚度，本发明在计算效率上得到了很大的提高，从原来的10个小时左右缩短至2-3分钟。

附图说明

图1为本发明的实现流程图；

图2为本发明实施例的直齿圆柱齿轮受力的示意图；

图3为本发明实施例的直齿圆柱齿实体齿基结构的示意图；

图4为本发明实施例的直齿圆柱齿只含有腹板的齿基结构的示意图；

图5为本发明实施例的直齿圆柱齿含有腹板和腹板孔的复杂齿基的示意图；

图6为本发明实施例的含复杂齿基的直齿圆柱齿轮有限元模型的示意图；

图7为本发明实施例的齿对活动线的示意图；

图8为本发明实施例的时变啮合刚度结果对比图。

具体实施方式

下面结合附图，并以表1所示的直齿圆柱齿轮结构参数为例，对本发明作进一步详细说明。

表1

参考图1，基于解析-有限元的复杂齿基直齿圆柱齿轮啮合刚度计算方法，包括如下步骤：

步骤1)基于解析法，分别计算直齿圆柱齿轮轮齿刚度，包括轮齿弯曲刚度值kb、剪切刚度值ks和径向压缩刚度值ka；

参考图2，给出了一个齿的示意图，将啮合力ft进行正交分解，得到轮齿上的轴向力fta和径向力ftb，其中，为接触点处的压力角，rb为基圆半径，cd段为过渡曲线部分，bc为渐开线曲线部分。

由于力对主动齿轮和从动齿轮的作用效果相同，故只取其中的一个齿轮进行分析，参考图2，径向力ftb会对轮齿产生弯曲和剪切作用，轴向力fta会对轮齿产生径向压缩作用，利用《材料力学》中的能量方法，计算轮齿相应的弯曲刚度值kb、剪切刚度值ks和径向压缩刚度值ka，具体计算公式表达式为：

其中，kb为轮齿的弯曲刚度，ks为轮齿的剪切刚度，ka为轮齿的径向压缩刚度，e为弹性模量，g为剪切模量，为分度圆压力角，rb为基圆半径，τc＝αc-(θb-invαc)，为接触点处的压力角，θb为基齿角的一半，为齿顶高系数，m为模数，γ和τ为积分范围，y1为过渡曲线上任一点处的水平坐标，y2为渐开线上任一点处的水平坐标，iy1为过渡曲线上任意位置处的截面惯性矩，iy2为渐开线上任意位置处的截面惯性矩，ay1为过渡曲线上任意位置处的截面面积，ay2为渐开线上任意位置处的截面面积。

步骤2)基于赫兹接触理论，计算轮齿的接触刚度值kh；

由于主动直齿圆柱齿轮与从动直齿圆柱齿轮在传动过程中彼此之间存在接触，因此存在接触刚度，采用hertz公式和能量法进行计算，其接触刚度的表达式为：

其中，kh为轮齿之间的接触刚度，e为弹性模量，b为切片的宽度，v为切片的泊松比。

步骤3)基于刚度串并联理论，计算单齿对啮合刚度值kt；

在计算单齿对的啮合刚度时，需要分别计算主动齿轮轮齿的弯曲刚度kb1、剪切刚度ks1、径向压缩刚度ka1和从动齿轮轮齿的弯曲刚度kb2、剪切刚度ks2、径向压缩刚度ka2以及轮齿之间的接触刚度kh，所以单齿对的啮合刚度kt是由主动齿轮轮齿的刚度、从动齿轮轮齿的刚度和轮齿之间的接触刚度kh串联构成，根据刚度的串并联关系，得到单齿对的啮合刚度值kt，其计算表达式为：

其中，kt为单齿对的啮合刚度，kb1为主动轮轮齿的弯曲刚度，ks1为主动轮轮齿的剪切刚度，ka1为主动轮轮齿的径向压缩刚度，kb2为从动轮轮齿的弯曲刚度，ks2为从动轮轮齿的剪切刚度，ka2为从动轮轮齿的径向压缩刚度，kh为轮齿之间的接触刚度。

步骤4)根据含复杂齿基的直齿圆柱齿轮几何参数建立含复杂齿基直齿圆柱齿轮的几何模型，并基于有限单元法，采用ansys软件中的solid185单元进行网格映射划分，形成含复杂齿基的直齿圆柱齿轮的有限元模型；

参考图3到图5，常见的齿轮基体主要有3种结构，图3表示实体齿基，图4表示只含有腹板的齿基，图5表示含有腹板和腹板孔的复杂齿基，也是本发明面向的情况。参考图6，,将直齿圆柱齿轮简化为一个圆柱体，圆柱体的高度等于直齿圆柱齿轮的齿宽，圆柱体的半径等于直齿圆柱齿轮的齿根圆半径，并根据直齿圆柱齿轮的参数，在圆柱体上完成腹板部分的建模和在腹板上完成减重孔的建模，从而形成含复杂齿基的直齿圆柱齿轮的几何模型，并使用solid185单元进行网格映射划分，其中设置单元的弹性模量为2.1×10¹¹pa，泊松比为0.3。

步骤5)根据齿轮啮合理论，建立刚性耦合区域；

步骤5.1)基于直齿圆柱齿轮的啮合理论，分别在主动齿轮和从动齿轮的啮合点位置处建立关键点，并采用mass21单元划分关键点，形成主质量点；

参考图6，根据直齿圆柱齿轮的啮合理论，计算出主动齿轮和从动齿轮的啮合点位置，并在位置处分别建立关键点，并使用ansys软件中的mass21单元划分关键点，形成主质量点。

步骤5.2)基于主质量点与复杂齿基上轮齿的对应关系，采用mpc184单元将主质量点与对应的齿基部分进行刚性耦合，形成刚性耦合区域；

参考图6，使用ansys软件中的mpc184单元，将步骤5.1)中建立的主质量点与步骤4)中建立的复杂齿基模型中对应齿基部分的节点进行连接，形成刚性耦合区域，从而模拟直齿圆柱齿轮的轮齿部分。

步骤6)根据齿轮啮合理论，判断主质量点对应的轮齿是否处于啮合状态，若处于啮合状态，则将处于啮合状态的轮齿对应的主质量点用ansys软件中的combin14单元连接，并将ansys软件中的combin14单元的刚度值设置为步骤3)中得到的kt值；

参考图7，根据啮合理论，首先判断步骤5.1)中建立的主质量点是否处于啮合状态，n1n2为理论活动线，b1b2为实际活动线，c1c2为建立的主质量点。n1c1和n2c2计算如下：

其中，mod函数为取模函数，z1为主动齿轮的参数，ω1t为主动齿轮的旋转角度，rb1为主动齿轮的主动轮中心孔半径。

若主质量点处于啮合状态，参考图6，采用ansys软件中的combin14单元将主动齿轮和从动齿轮上对应的主质量点进行连接，并将ansys软件中的combin14单元的刚度值设置为步骤3)中得到的kt值；

步骤7)基于静态分析，计算直齿圆柱齿轮的时变啮合刚度值k。

参考图6，对模型施加约束，首先在主动齿轮中心孔处施加扭矩t；其次在主动齿轮中心建立主质量点，并将该主质量点与主动齿轮内孔进行自由度耦合，同时约束该主质量点除绕自身轴线转动的其余自由度；然后采用相同的办法在从动齿轮中心建立主质量点，并将该主质量点与从动齿轮的内孔进行自由度耦合，同时约束该主质量点的全部自由度；最后采用ansys的静态分析对模型进行计算，得到模型的变形结果，通过提取主动齿轮中心的主质量点绕轴线的转动变形量δθ，得到直齿圆柱齿轮的时变啮合刚度值k，其计算表达式为：

其中，t为施加在主动轮中心孔处的扭矩，rb1为主动轮中心孔半径，δθ为主动轮中心孔的变形。

步骤8)判断是否完成一个啮合周期的刚度计算，若没有完成一个啮合周期的刚度计算，则直齿圆柱齿轮旋转一定角度，重复步骤1)至步骤7)，若完成了一个啮合周期的刚度计算，则停止计算，得到含复杂齿基直齿圆柱齿轮的时变啮合刚度。

基于上述步骤，结合本发明实施例的具体参数进行计算，参考图8，图中实线表示按照本发明所述方法计算得到的直齿圆柱齿轮时变啮合刚度曲线，虚线表示只采用有限元法进行模拟得到的直齿圆柱齿轮时变啮合刚度曲线，得到本发明实施例的时变啮合刚度的最大相对误差为2.03％，并且计算时间由传统方法的10h左右缩短至2-3min，相比只采用有限元法模拟的现有技术，计算效率得到大幅提高。