基于ANSYS workbench齿轮啮合刚度计算及动力学仿真*

何育民，郝安帮

(西安建筑科技大学 机电工程学院，西安 710055)

齿轮是机械装置中的重要零部件，广泛应用于汽车、航空等领域.啮合刚度是齿轮系统振动的内部激励源之一，正确地求解啮合刚度以及齿轮啮合动力学研究对工程实际具有重要意义.

王旭等[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

图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 轮齿接触设置

图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中的MassCoefficient和StiffnessCoefficient.在大多数实际的结构问题中，α阻尼可以被忽略，即α=0，此时β值由式(2)确定，即

β=2εi/ωi

(2)

式中：εi为阻尼比，根据实验值一般为0.03～0.17，这里取0.03；ωi为第i阶固有频率，通过有限元模态分析得到ω1=646.33Hz.图5是齿轮啮合的前六阶模态阵型图.从图5中可看出一阶模态阵型是两齿轮相对反向旋转，二阶模态阵型是两齿轮相对同向旋转，其他阶数模态阵型则不是齿轮副的轴向转动.而实际齿轮啮合过程中，主从动轮是相对反向旋转的，且主要是此模态下的振动，故取一阶模态频率，所以β=4.64×10-5.NumericalDamping这个选项控制着一个结构较高频率所引起的数值噪声，通常这些高频模态的贡献是不精确的，为了避免这种现象，给定一些数值阻尼是可取的，对于瞬态结构分析，数值阻尼设为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′thBMBRBBcosβ

(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) 对不同轮齿摩擦系数的直齿圆柱齿轮进行动力学仿真，发现轮齿摩擦系数的增大会导致齿轮角加速度和齿轮旋转副竖直与水平约束力的最大幅值增大，轮齿摩擦系数对齿轮副的轴向约束力的影响比较显著.