摘要：本文使用高性能并行计算机和ANSYS软件，计算一个真实齿轮的三维齿根应力和变形，为齿根应力和轮齿变形提供了一个计算标准，为齿轮动态、优化和可靠性设计和CAE提供新的基础，也测试了ANSYS软件的有限元并行计算能力。计算是在曙光4000A超级计算机上进行，计算模型的自由度为659万，工况5个，计算时间为2个小时，使用了4个CPU。
关键词：齿根应力，轮齿变形，超级计算机，并行计算，计算标准
1 前言
齿轮齿根应力和轮齿变形的计算是齿轮设计和研究的基本问题。过去由于计算机的限制，许多研究和分析人员计算齿根应力和轮齿变形时所用的模型为二维或三维等效轮齿模型，计算模型存在误差，计算模型和结果的可靠性低，因此计算结果很难比较，特别是与试验结果不好比较，造成很多重复试验。所以有必要建立一个三维整轮齿根应力和轮齿变形计算的标准，给齿轮设计和研究人员和CAE研究人员参考。使用高性能计算，现在能够对齿轮三维整轮模型进行仿真分析，精确计算齿轮齿根应力和轮齿变形。本文对一个采煤机齿轮进行了计算，计算是在曙光4000A超级计算机上进行。
2 齿轮齿根应力和轮齿变形的计算
表1为计算齿轮和加工刀具的参数，刀具为凸头滚刀，齿顶高系数为1，齿轮配合的轴直径（齿轮内圆直径）为90毫米。表 1 计算齿根应力齿轮参数
 2.1 齿轮三维计算模型
载荷均匀或线性分布作用在齿顶，作用载荷大小设定为F = 100N/mm，不同工况总载荷大小一样，齿轮内圆固定约束。齿轮弹性模量E = 2.06e5 N/mm2，泊松比ν=0.3。
2.1.1 网格模型
图1是计算齿轮实体几何模型，由于齿轮齿根有应力集中，应该选择三维二次等参单元，ANSYS软件中使用单元SOLID95表示，齿轮整轮网格如图2所示。整轮网格划分分三部分，载荷作用轮齿扇、相邻轮齿扇（左、右）和其余轮齿扇，每部分分轮齿和轮缘（轮幅）两部分。载荷作用轮齿扇的轮齿部分（图3）由于齿根网格密度高，因此其它部分网格密度高也相应提高，轮缘部分周向网格密度最高（相对其它轮缘部分）。由于应力主要集中在载荷作用轮齿扇局部范围内，因此与它间隔轮齿扇可以使用较稀疏的网格，相邻轮齿扇的网格是载荷作用轮齿扇和间隔轮齿扇之间的过渡网格。其余轮齿扇中的各扇网格一样，可以复制。所有轮缘部分径向网格密度一样，整轮网格齿向密度一样。载荷分布作用在齿顶的节点上，载荷作用轮齿的网格和均匀分布载荷如图3所示。单元网格三边长度比例有一定要求。 
图1 齿轮实体三维几何模型 图２ 齿轮三维整轮网格模型
图3 齿轮轮齿三维网格模型
2.1.2 计算时间和并行性能
齿轮计算模型有518940单元，2196347节点，共6589041自由度。计算工况为5个。计算时间在表2中，求解方法使用ANSYS的预条件迭代求解器。不同CPU数求解时间如图4所示，2个CPU的加速比为2.17，4个CPU的加速比为３.２２。ANSYS软件共享内存并行求解效率很高，不比分布式并行求解差，4个CPU的加速比可达３以上。表2 多工况三维整轮齿根应力有限元模型在曙光4000A上的计算时间
 
图4 计算求解时间（ＣＰＵ和墙钟ＷＡＬＬ时间）
2.2 齿轮齿根应力
不同工况和偏载方式的最大齿根应力计算结果在表3中，σf是受拉侧最大齿根拉应力（即最大主应力σ1），σC是受压侧最大齿根压应力（即最小主应力σ3值），表4是作用在齿顶的5种载荷工况。图5是均匀分布载荷受拉侧齿根应力齿向分布，σI =σ1 - σ3= 2ｔmax，为应力强度等于2倍最大切应力ｔmax，σE 为等效应力。对于平面应力或平面应变，理论上σ1=σI ，实际上由于计算误差不是绝对相等[2]。对于平面应力，理论上等效应力σE =σ1，而平面应变，理论上等效应力： 。根据图5计算结果所知，齿轮在中部大部分（60%）都工作在平面应变状态，而在最端部工作在平面应力状态。最大应力是在应变状态向平面应力过渡处，离端部距离大约是齿长的五分之一（23%）处，端部应力最小，只有最大应力的65%，中部应力比最大应力略小（小1.73%）。工况2，三倍载荷均匀分布作用在齿轮齿面中间三分之一处，但齿根应力仅增加2.1倍。工况3，三倍载荷均匀分布作用在齿轮齿面一端三分之一处，齿根应力增加了2.65倍。工况4，一端二倍载荷线性分布作用在齿轮齿面上，齿根应力仅增加1.63倍。工况5，一端三倍载荷线性分布作用在齿轮齿面上，齿根应力增加2.23倍。
图5 均匀分布载荷齿根应力齿向分布
σ1为最大主应力，σI =σ1-σ3 σE 为等效应力
表3 不同载荷工况的齿根应力及偏载系数
 表4 5种载荷工况
 2.3 齿轮轮齿变形
齿轮变形包括轮齿接触变形、挠曲变形和本体变形。轮齿变形文献[3]已经研究，但是仅仅是二维等效模型。本文研究计算了轮齿挠曲变形和本体变形三维模型的影响，同时也计算了不同载荷工况和偏载荷下的齿轮轮齿挠曲变形和本体变形。图6是工况1均匀分布载荷下的齿轮变形，u1是载荷作用轮齿对称点载荷作用线方向的变形，u 2是与载荷作用轮齿载荷作用线方向相临的轮齿对称点的变形。因此u 2近似表示齿轮本体变形，u 1- u 2可定义为轮齿挠曲变形。计算结果表明，均匀分布载荷齿轮变形沿齿向分布有点变化，u1变化与齿根应力一致，但幅度很小，u 1最大为5.461微米，最小5.321微米。u 2变化与u 1相反，，最大为1.460微米，最小1.444微米。因此齿轮轮齿变形挠曲为3.939微米，齿轮本体变形1.452微米。
图6 均匀分布载荷齿轮变形齿向分布
u1为载荷作用轮齿对称点变形，u2 为相临轮齿对称点变形
2.4 计算模型误差分析
由图2计算模型的离散误差已经很小，由图5计算结果，根据最大主应力和应力强度之差，齿根应力的计算误差估计在3‰。ANSYS软件计算的齿根应力局部离散误差仅为0.35‰，因此计算模型的可靠性很高，计算误差肯定小于1%。另一方面，齿轮计算模型的约束假定与轴配合的内圆为固定约束，实际上约束有多种，因而不同约束对齿根应力和齿轮变形肯定有影响。这个影响对齿根应力和轮齿变形肯定很小，对于轮缘厚度较大的，可以忽略（如本文），但是是否小于1%，还需要进一步研究。齿轮接触载荷是假定为均匀分布或者线性分布作用在齿面上，实际上接触载荷不是按照假定的理想分布，因此建立一对齿轮三维接触仿真分析模型是重要的课题。
齿轮变形计算模型的离散误差比齿根应力小，计算误差肯定小于1%，但齿轮变形齿向分布的变化可能要大于1%。本文均匀载荷齿轮本体变形齿向分布的变化是1.1%，总变形变化2.6%（图6），因此计算轮齿变形可能的误差有1%左右。齿轮变形主要受弹性模量影响，弹性模量的误差和随机性，可能使得齿轮变形的误差达到1-5%，需要进一步研究。另一方面，齿轮变形（特别是本体变形）受齿轮本体结构影响很大，本体结构可能使得齿轮变形的误差达10%以上[4]，但对轮齿变形影响较小，可能也会达到1%，因此也需要进一步研究。同样不同载荷分布齿轮变形沿齿向是变化的，需要建立一对齿轮三维接触仿真分析模型。
3 结论
通过建立三维整体齿轮仿真分析模型，使用高性能计算机，可以精确计算齿根应力和齿轮变形，计算模型的计算误差肯定小于1%。 不同载荷分布下齿根应力和轮齿变形的分布能够精确计算，齿根应力和轮齿变形的一个新的计算标准已经建立，三维模型对齿根应力和轮齿变形的影响已经显示，齿根应力和轮齿变形计算得到了进一步发展，齿轮动态设计、优化设计、可靠性设计和CAE有了新的基础。
ANSYS软件的并行计算得到应用和检验，超级计算得到应用，659万个自由度成功在曙光4000A上有效运行，并行计算的优势得到发挥，新的ANSYS并行计算基准产生。
[参考文献]
[1] 吴继泽，王统编著﹒齿根过渡曲线和齿根应力﹒国防工业出版社，1989。
[2] 杨生华，汪勤悫﹒新的有限元法在齿根应力计算中的应用和趋势﹒“力学2000年”学术大会论文集，2000.8: 507-509。
[3] 杨生华﹒齿轮接触有限元分析. 计算力学学报，2003，20（2）：189-194。
[4] 杨生华.子结构、子模型方法在齿轮应力和变形计算中的应用. 工程与科学中的计算力学， 北京:北京大学出版社，2001: 701-707。
[5] 杨生华﹒高性能计算和CAE﹒第二界中国CAE工程分析技术年会（CCATAX，2006）论文集，2006.7。(end)
|
