齿轮是机械装置中的重要零部件，广泛应用于汽车、航空等领域.啮合刚度是齿轮系统振动的内部激励源之一，正确地求解啮合刚度以及齿轮啮合动力学研究对工程实际具有重要意义.
王旭等[1]用材料力学的方法计算了正常和含裂纹齿轮的啮合刚度，并且对有、无裂纹齿轮副进行了动力学特性分析；万志国等[2]考虑了齿轮基体变形与齿根圆、基圆不重合的情况，进一步修正了势能法求解啮合刚度的公式；Zhan等[3]利用有限元软件workbench用准静态方法求解了啮合刚度，并研究了齿顶圆角及主从动轮轴不平行时对啮合刚度的影响；毕凤荣等[4]在MSC-ADMAS平台上对渐开线齿轮副动态啮合力进行了计算机仿真，并研究了转速对动态啮合力的影响；陈占锋等[5]在Pro/Engineer软件中建立了两级行星齿轮三维模型，并导入ADMAS软件中进行动力学仿真，验证了行星齿轮系统的传动规律，为其强度、刚度和疲劳分析提供了依据；吴勇军等[6]提出了一种可综合考虑齿轮副连续弹性啮合过程中多种影响因素的接触有限元分析方法.
啮合刚度是齿轮动力学分析中的重要参数，而齿轮副在不同的工况下，实际重合度是不相同的，齿轮动力学数值仿真时需要啮合刚度有明确的数学表达式，因此多数学者对啮合刚度进行了曲线拟合.引用最多的是文献[7]的方法，这些公式中都包含重合度的因素，若将理论重合度下的啮合刚度公式代入齿轮动力学方程，将会因为理论重合度与实际重合度不同，导致动力学分析出现偏差.
本文在三维建模软件Solidworks中利用渐开线参数方程作为曲线驱动方程，建立了直齿圆柱齿轮副装配体模型，然后导入ANSYS workbench瞬态动力学模块，用准静态方法求出啮合刚度，并与ISO6336-1-2006啮合刚度计算值比较，验证了模型建立的精确性以及有限元设置的合理性，同时也验证了实际工况下因为重合度的改变，理论重合度计算的啮合刚度与实际工况下啮合刚度不同步的现象.
1 齿轮副模型的建立及有限元设置
1.1 齿轮副模型建立与网格划分
本文所建立的直齿圆柱齿轮相关参数如表1所示，根据表1信息在Solidworks里建立直齿圆柱齿轮副三维模型.图1为建立的单个齿轮模型图，其中，Rb、Rf、r、p分别为齿轮基圆半径、齿根圆半径、过渡圆弧半径及分度圆齿距.过渡圆角半径r=0.38m，m为模数.齿轮啮合部分为轮齿，为了减少网格数量，齿轮轴孔半径取20 mm.
表1 直齿圆柱齿轮参数
Tab.1 Spur gear parameters
齿数模数mm齿厚mm压力角(°)齿宽mm变位系数中心距mm密度(kg·m-3)杨氏模量MPa泊松比重合度3029.4201006078502000000.31.6535
图1 齿轮啮合三维模型
Fig.1 3D model for gear meshing
齿廓曲线选择方程式驱动的曲线命令进行创建，输入渐开线参数方程式，剪去多余的线条，对单个齿廓进行圆周阵列，然后轴向拉伸就得到单个齿轮，最后按照齿轮副中心距进行装配.
将齿轮副装配体模型通过Solidworks与ANSYS workbench的接口直接导入ANSYS workbench的瞬态动力学模块，由于每对轮齿啮合有相同的规律，为了缩短有限元求解时间，故只保留8对轮齿.图2是主从动轮施加边界条件的过程，左边是主动轮，右边为从动轮，主从动轮轴孔处添加旋转副约束.主动轮旋转副施加10 r/min恒定转速，从动轮旋转副施加200 N·m转矩.由于参与啮合的主要是轮齿，所以轮齿网格加密，齿轮基体网格可以稍微稀疏一些，图3是齿轮副整体网格.旋转副使用MPC184单元，接触区使用CONTA174和TARGE170单元，因为二次单元会导致等效节点接触力在角节点和边节点之间震荡，所以齿轮副主体使用solid185单元，共得到178 220个单元，197 243个节点.
1.2 轮齿接触设置
图4是有限元软件ANSYS workbench瞬态动力学接触设置与分析设置的界面.图4a为接触设置界面.啮合的轮齿定义为面
面接触，而每个轮齿只有一个面参与啮合，故有8个面面接触对.法向刚度因子越小越容易收敛，但是会引入大的穿透值，这将导致求解结果不精确，法向刚度因子的确定方法见文献[8]，通过不断调整最终设为1，使用扩展拉格朗日算法可以减少求解过程中对刚度的敏感性.在接触初始阶段，接触与目标单元积分点之间有微小间隙存在，导致接触不会被检测到进而造成接触体之间的刚体运动.但对于非线性接触来说，尽管忽略掉起始间隙，间隙仍会在加载过程中出现，Adjust to Touch的作用是将即将发生的接触限制在检测的范围内，如果接触前夕有间隙和穿透存在，将忽略间隙和穿透直接接触，其他设置保持默认.
图2 约束设置与载荷施加
Fig.2 Constraint setting and load implementation
图3 整体网格
Fig.3 Whole meshing
图4 接触与分析设置
Fig.4 Settings for contact and analysis
1.3 分析设置
图4b是分析设置界面.将Time Integration设置为off，则workbench瞬态动力学求解中不包括惯性效应和热容，不计算速度结果，因此由速度产生的弹簧阻尼力将为0，这是一个准静态过程.阻尼控制有两种，这里选择瑞利阻尼，其表达式为
C=αM+βK
(1)
式中：C为阻尼矩阵；M为质量矩阵；K为刚度矩阵；α、β为阻尼常数，分别对应图4b中的Mass Coefficient和Stiffness Coefficient.在大多数实际的结构问题中，α阻尼可以被忽略，即α=0，此时β值由式(2)确定，即
β=2εi/ωi
(2)
式中：εi为阻尼比，根据实验值一般为0.03～0.17，这里取0.03；ωi为第i阶固有频率，通过有限元模态分析得到ω1=646.33 Hz.图5是齿轮啮合的前六阶模态阵型图.从图5中可看出一阶模态阵型是两齿轮相对反向旋转，二阶模态阵型是两齿轮相对同向旋转，其他阶数模态阵型则不是齿轮副的轴向转动.而实际齿轮啮合过程中，主从动轮是相对反向旋转的，且主要是此模态下的振动，故取一阶模态频率，所以β=4.64×10-5.Numerical Damping这个选项控制着一个结构较高频率所引起的数值噪声，通常这些高频模态的贡献是不精确的，为了避免这种现象，给定一些数值阻尼是可取的，对于瞬态结构分析，数值阻尼设为0.1.
2 啮合刚度计算
图6是不计轮齿摩擦力的齿轮副动力学模型图.图6中主动轮受到驱动力矩与啮合力而逆时针方向旋转；从动轮在啮合力与负载转矩共同作用下顺时针旋转.动力学方程为
(3)
(4)
(5)
式中：φp、φg为主、从动轮的转角；Ip、Ig为主、从动轮的轴向转动惯量；rpb、rgb为主、从动轮的基圆半径；cm(t)为啮合阻尼；km(t)为齿轮时变啮合刚度；x为动态传递误差；e(t)为静态传递误差；Tp、Tg分别为外界对主从动轮的驱动力矩.
瞬态动力学模块在准静态环境下，不包括齿轮副惯性力以及阻尼力，由渐开线建立的齿轮副模型标准安装，可认为静态传递误差e(t)=0，由式(3)～(5)可得啮合刚度表达式为
(6)
图5 前六阶模态振形
Fig.5 Vibration shapes of first six order modals
通过仿真求出主从动轮转角φp和φg，代入式(6)就可求出综合啮合刚度.
图6 齿轮啮合动力学模型
Fig.6 Dynamics model for gear meshing
仿真完成后，提取主从动轮的转角，得到动态传递误差，计算出啮合刚度.图7是有限元计算的啮合刚度与文献[7]计算的啮合刚度对比图.从图7中可以看出，理论重合度计算的综合啮合刚度与实际载荷有限元计算的啮合刚度出现了不同步的现象，而大多数齿轮动态响应是在理论重合度计算的综合啮合刚度上进行的.轮齿啮合具有周期性，为了直观地体现这种规律，引入无量纲啮合时间t/tc，t为绝对时间，tc为轮齿延啮合线移动一个分度圆法向齿距的时间.
图7 啮合刚度对比图
Fig.7 Meshing stiffness contrast diagram
ISO 6336-1-2006标准中规定，对于重合度ε≥1.2，螺旋角β≤30°的斜齿轮，其单齿啮合刚度最大值为
B′=B′th BM BR BBcos β
(7)
式中：B′th为理论单齿刚度；BM为理论修正系数；BR为轮胚结构系数；BB为基本齿廓系数.
ISO6336-1-2006方法求得的单齿最大啮合刚度和平均啮合刚度分别为130 931.5和195 104.3 N/mm，本文计算的单齿最大啮合刚度和平均啮合刚度分别为137 314.4和205 054.6 N/mm，相对误差分别是4.9%和5.1%，这验证了所建立的齿轮副模型的精确性以及有限元设置的合理性.
3 轮齿啮合动力学仿真
在前面建立的有限元模型基础上，考虑齿轮副的惯性效应和阻尼效应，将图4b中time integration调整为on，完成直齿圆柱齿轮副动力学仿真.图8、9分别是从动轮角速度、角加速度仿真图；图10、11分别是主动轮旋转副处的水平与竖直约束力仿真图.根据啮合轮齿接触面压力的出现与消失可以划分单双齿啮合区，结合图8～11可以看出，齿轮开始进入和退出双齿啮合时振动加剧，轮齿的交替啮合产生冲击力，单齿啮合比双齿啮合振动强烈.
图8 从动轮角速度
Fig.8 Angular velocity of driven wheel
图9 从动轮角加速度
Fig.9 Angular acceleration of driven wheel
图10 主动轮水平约束力
Fig.10 Horizontal constraint force of driving wheel
轮齿摩擦系数不仅影响齿形的相对拉尖程度[9]，而且也影响着齿轮的啮合.在前面模型的基础上，进一步考虑啮合轮齿之间的摩擦力，分别将摩擦系数设为0.02、0.03进行动力学仿真，得到从动轮的加速度以及主动轮旋转副处的约束力数据.为了直观地看出变化规律，分别提取啮合稳定后的100组数据在MATLAB中进行快速傅里叶变换，得到不同轮齿摩擦系数下从动轮的角加速度频域图(如图12所示)及主动轮旋转副处的水平、竖直以及轴向约束力频域图(如图13～15所示).
图11 主动轮竖直约束力
Fig.11 Vertical constraint force of driving wheel
图12 加速度幅值频域信号
Fig.12 Acceleration amplitude signal in frequency domain
图13 水平约束力幅值频域信号
Fig.13 Horizontal constraint force amplitude signal in frequency domain
图14 竖直约束力幅值频域信号
Fig.14 Vertical constraint force amplitude signal in frequency domain
图15 轴向约束力幅值频域信号
Fig.15 Axial constraint force amplitude signal in frequency domain
仿真中主动轮施加的转速为10 r/min，转频约为0.17 Hz，啮合频率为5 Hz.由图12可知，不同摩擦系数的从动轮角加速度最大幅值均出现在50 Hz的频率处(是啮合频率的10倍)，且随着摩擦系数的增大而增大；由图13、14可知，约束力最大幅值出现在1 Hz的频率处，是转频的6倍频，且最大幅值随轮齿摩擦系数的增大而增大；由图15可知，虽然直齿轮副啮合时轴向约束力很小，但是轮齿摩擦系数对直齿圆柱齿轮副轴向约束力的影响比较明显.
4 结 论
通过上述分析可以得出以下结论：
1) 在三维建模软件Solidworks中利用渐开线参数方程建立了直齿圆柱齿轮副装配体模型，并导入ANSYS workbench的瞬态动力学模块中，用准静态的方法求出啮合刚度，与ISO6336-1-2006计算的啮合刚度比较吻合.实际工况下综合啮合刚度与理论重合度计算的综合啮合刚度产生了不同步的现象.
2) 通过直齿圆柱齿轮动力学仿真可以发现，轮齿在开始进入和退出双齿啮合时振动加剧，轮齿的交替啮合会产生冲击力，单齿啮合比双齿啮合振动强烈.
3) 对不同轮齿摩擦系数的直齿圆柱齿轮进行动力学仿真，发现轮齿摩擦系数的增大会导致齿轮角加速度和齿轮旋转副竖直与水平约束力的最大幅值增大，轮齿摩擦系数对齿轮副的轴向约束力的影响比较显著.
[1]王旭，伍星，肖正明，等.含裂纹故障的齿轮系统动力学特性研究及其故障特征分析 [J].振动与冲击，2017，36(9)：74-79.
(WANG Xu，WU Xing，XIAO Zheng-ming，et al.Dynamic characteristics of a gear system with crack fault and its fault feature analysis [J].Journal of Vibration and Shock，2017，36(9)：74-79.)
[2]万志国，訾艳阳，曹宏瑞，等.时变啮合刚度算法修正与齿根裂纹动力学建模 [J].机械工程学报，2013，49(11)：153-160.
(WAN Zhi-guo，ZI Yan-yang，CAO Hong-rui，et al.Time-varying mesh stiffness algorithm correction and tooth crack dynamic modeling [J].Journal of Mechanical Engineering，2013，49(11)：153-160.)
[3]Zhan J X，Fard M，Jazar R.A CAD-FEM-QSA integ-ration technique for determing the time-varying me-shing stiffness of gear pairs [J].Measurement，2017，100：139-149.
[4]毕凤荣，崔新涛.渐开线齿轮动态啮合力计算机仿真 [J].天津大学学报(自然科学与工程技术版)，2005，38(11)：991-995.
(BI Feng-rong，CUI Xin-tao.Computer simulation fordynamic meching force of involute gears [J].Journal of Tianjin University (Science and Technology)，2005，38(11)：991-995.)
[5]陈占锋，程刚，陈曦晖，等.两级行星齿轮动力学仿真研究 [J].现代制造工程，2015(12)：70-73.
(CHEN Zhan-feng，CHENG Gang，CHEN Xi-hui，et al.Research on dynamics simulation of two-stage planetary gear [J].Modern Manufacturing Engineering，2015(12)：70-73.)
[6]吴勇军，梁跃.齿轮副动态啮合特性的接触有限元分析 [J].振动与冲击，2012，31(19)：61-67.
(WU Yong-jun，LIANG Yue.Dynamic meshing characteristics of gear pair using contact finite element method [J].Journal of Vibration and Shock，2012，31(19)：61-67.)
[7]张建云，丘大谋.一种求解直齿圆柱齿轮啮合刚度的方法 [J].西安建筑科技大学学报，1996，28(2)：134-137.
(ZHANG Jian-yun，QIU Da-mou.Method for calcu-lating spur gear mesh stiffness [J].Journal of Xi’an University of Architecture and Technology，1996，28(2)：134-137.)
[8]白恩军，谢里阳，佟安时，等.考虑齿轮轴变形的斜齿轮接触分析 [J].兵工学报，2015，36(10)：1975-1981.
(BAI En-jun，XIE Li-yang，TONG An-shi，et al.The contact analysis of helical gear in considering gear shaft deformation [J].Acta Armamentari，2015，36(10)：1975-1981.)
[9]朱小星，王宝雨，付晓斌.热轧齿轮齿形拉尖的影响因素 [J].沈阳工业大学学报，2016，38(4)：410-415.
(ZHU Xiao-xing，WANG Bao-yu，FU Xiao-bin.Influencing factors for tooth tip pulling of hot roll forming gear [J].Journal of Shenyang University of Techno-logy，2016，38(4)：410-415.)