林支慨，李瑞阁，燕 楠

(1.上海理工大学机械工程学院，上海200093； 2.上海理工大学光电信息与计算机工程学院，上海200093)

0 引言

接触情况在机械工程中是广泛存在的，尤其是在运动部件中，如轴承滚动体与内外圈的接触、轨道车辆与钢轨间的接触、轮胎与地面的接触。其中齿轮接触是接触情况中典型的问题，同时对于很多学者来说也是一大难点。许多学者为了解决一些齿轮问题花费了很多时间去研究，并且取得了一些成果。付昆昆等[1]用石川法和改进石川法来计算齿轮多齿时变啮合的刚度数值，并与有限元法进行了对比，结果表明，改进石川法更能反映真实情况，结果更加精确。MAO K[2]使用微几何修正法降低了齿面疲劳磨损程度，使得有限元计算的接触应力更加精确。李学志等[3]为了解决斜齿轮在啮合更替过程中受到载荷突变造成齿面点蚀和冲击的问题，对斜齿轮接触问题和修形问题进行了分析研究。卓耀彬等[4]为了使准双曲面齿轮传动系统更适合静态齿面接触分析，运用自己的方法重新建立了准双曲面齿轮传动系统的有限元模型。黄学文[5]为了能够精确地得出斜齿轮副接触时的状态和应力大小的一般规律，用B样条曲面在ANSYS中建立了精确的斜齿轮接触齿面实体模型，并用有限元法对模型进行了接触分析。

齿轮在实际应用中失效的主要形式是齿面损坏，所以在齿轮设计过程中对齿轮进行精确地分析是很重要的一个阶段。这些接触问题不仅关系到所设计出的产品寿命及功能是否满足需求，甚至有可能会危及生命财产安全[6-7]。随着科技的发展，齿轮的精度要求也更加严格。除了像设计方面的内在因素外，外在因素也是影响齿轮品质的重要原因，如齿轮热胀时啮合会产生几何干涉，而长时间在热胀区工作时会造成齿面胶合失效等问题，这些会直接导致齿轮失效，进而造成事故。目前，国内在接触强度理论和弯曲理论方面的研究也有很多，如李广兴[8]为了解决某种汽车变速器二挡时齿轮副发出噪声的问题，进行了噪声分析仿真和一系列试验研究。为了预测出齿轮和轴承的寿命，他们对齿轮激励源和NVH问题进行了研究和分析，最后提出了新的修改方案，并解决了变速器噪声问题。高准等[9]为了能够精确计算齿轮在特定载荷下的变形情况，对基于Hertz接触理论模型的行星传动齿轮接触问题进行了分析研究。

目前，有限元法是解决接触分析问题中最常用的方法，在工程分析方面有限元法也是最经典的方法之一。它适用于多种领域，如弹塑性力学、热传导、流体力学等领域。虽然有限元的概念早就被学者们提出，但是由于当时计算机还没有能够被人们应用，有限元法也就没有得到人们的重视。所以有限元的应用是离不开计算机发展的，因为计算机的不断发展和应用使得有限元法在工程领域的应用成为可能，也只有在计算机的基础上才会体现出有限元法的优势。有限元的核心思想是结构离散化，是将复杂的模型用有限个相同规则的单元进行表示，目前应用的单元只有四边形单元和六边形单元。对于非结构化有限元网格来说，四边形和六边形有限元网格由于在计算精度、划分网格数量等方面具有许多优势，已经成为有限元分析中首选网格[8]。

很多文献对于接触分析其实并不是很精确，完全属于理想状态，并没有考虑一些因素对模型接触分析的影响，导致分析的结果不精确。本文在常规有限元分析的基础上考虑表面切向应力对接触应力分布的影响，结果表明，表面切向应力不仅影响模型的接触应力分布，还会影响最大接触应力值。本文还比较了四边形与六边形有限元网格在分析时的效率问题，在相同的网格单元大小的前提下比较接触应力的收敛速度。结果表明，在相同大小的网格单元条件下，四边形网格的收敛速度要比六边形快，说明四边形网格的分析效率相对更高效一些。这些都为以后的分析工作提供了参考价值，提高了分析模型的精确性和效率。

1 有限元接触理论

弹性接触问题最早是由赫兹提出并进行研究的，他提出了经典的赫兹接触理论，但是该理论具有一定的局限性，在工程实际应用中只能对形状简单、接触面相对简单规整的问题进行分析和计算。因为这些条件的限制，赫兹理论主要被应用在齿轮接触分析、轴承滚动体接触分析、轨道车辆与钢轨的接触分析等领域。接触理论包括Hertz解析算法和数值解法。

赫兹问题的解析解可通过赫兹公式计算得到。根据经典赫兹理论，两个圆柱的最大接触应力的理论公式为

(1)

式中E1、E2——弹性模量

μ1、μ2——泊松比

b——厚度

ρ1、ρ2——半径

F——外力

接触面半宽

(2)

许用接触应力

(3)

式中σHlim——齿轮的接触疲劳极限

SH——齿面接触疲劳安全系数

Von Mises应力是基于剪切应变能的等效应力。它用等值线来表示模型内部的应力分布情况，这样可以清晰地表达出分析结果在整个模型中的变化，从而可以知道模型中最需要优化和改进的地方。求范式应力的理论公式为

σeq=

(4)

式中σ1，σ2，σ3——第一、二、三主应力

由弹性理论可以得到体积力与内力的关系方程

(5)

表面力与体内的应力存在的边界条件

(6)

式中FSx、FSy、FSz——表面力在x、y、z方向的分量

l、m、n——表面外法线方向的余弦

2 有限元接触分析

2.1表面切向应力对接触分析的影响

当接触变形条件和载荷条件确定时，都只会考虑表面的法向集中载荷，而忽略表面切向应力的影响。为了能够更加精确地分析出接触应力的分布，并且考虑零件表面疲劳的耐久性，本文考虑了表面切向应力。在有些应用中是可以不考虑表面切向力影响的，特别是在滚动轴承的应用，因为有充分的润滑，所以在运动过程中都会忽略它们之间的滑动。这就表明，切向应力与法向应力相比是非常小的，所以切向应力在这种情况是可以忽略的。

有些情况是不能忽略的，如在测某零件表面疲劳的耐久性时，表面切应力就不能被忽略。因为表面切应力对于特定情况下的零件表面耐久性的影响是很重要的，所用到的计算方法也是很复杂的。

在考虑表面切向应力的情况下，对于两接触物体在接触部位的von Mises应力的位置与法向应力和表面切应力的比值是有很大关系的。从图1中可以看出不同应力出现的深度，随着表面切应力与法向应力比值的增加，最大von Mises应力更接近于接触表面，当τ/σ=0.3时，最大von Mises应力会出现在接触表面上[10]。

图1 不同的表面切应力(τ/σ)和材料应力 (σVM/σ)随深度的变化Fig.1 Variation of surface shear stresses(τ/σ)and material stress(σVM/σ)with depth

2.2有限元接触分析实例

两个相同的1/4圆柱接触的受力分析如图2所示。圆柱半径0.5 mm，厚度8 mm。材料的弹性模量E=3.5×106MPa，泊松比μ=0.3。载荷位置为上半圆的上表面，大小为100 N。约束位置为下半圆的下表面。

图2 1/4圆柱接触Fig.2 Quarter-cylinder contact

对于六面体网格单元，单元大小为0.05 mm的1/4圆柱接触分析结果如图3～4所示，不同单元大小

图3 六面体网格单元1/4圆柱接触分析Fig.3 Contact analysis of a quarter cylinder for hexahedral mesh element

下的分析结果如表1所示。从分析结果可以看出，最大接触应力在模型接触表面上。

图4 六面体网格单元接触分析局部放大Fig.4 Local magnification of contact analysis for hexahedral mesh element

表1 六面体网格单元分析结果

对于四面体网格单元，单元大小为0.05 mm的1/4圆柱接触分析结果如图5～6所示，不同单元大小下的分析结果如表2所示。从分析结果可以看出，最大接触应力的分布与六面体的有所不同，接触表面和接触次表层同时出现了最大接触应力。

图5 四面体网格单元1/4圆柱接触分析Fig.5 Contact analysis of a quarter cylinder for tetrahedral mesh element

图6 四面体网格单元接触分析局部放大Fig.6 Local magnification of contact analysis for tetrahedral mesh element

表2 四面体网格单元分析结果

从表1～2中可以看出，在单元大小一样的条件下，四面体单元分析的结果要比六面体单元的分析结果更精确，从收敛曲线(图7)可以看出，四面体网格分析效率更高。为了使分析结果更加精确，可以使用局部画网格的方法。在知道接触应力最大值和位置情况下，可以对应力集中的地方进行局部细化，然后再进行分析，如图8～10所示。

异戊烯醇聚氧乙烯醚(TPEG)：分子量1200；丙烯酸(AA)；醋酸乙烯酯(VAc)；马来酸二乙酯(MADE)；丙烯酸羟丙酯(HPA)；烯丙基磺酸钠(ALS)；过氧化二苯甲酰(BPO)、偶氮二异丁腈(AIBN)、过硫酸铵(APS)。

图7 不同网格大小的最大接触应力Fig.7 Maximum contact stress of different mesh sizes

图8 四面体网格划分模型Fig.8 Tetrahedral meshing model

图9 1/4圆柱接触应力分布Fig.9 Contact stress distribution of a quarter cylinder

从分析结果可以看出，最大接触应力出现在模型接触表面的次表层，这是表面切应力对接触应力的影响，最大接触应力为1 931.7 MPa。根据赫兹理论的解析解公式可以计算得最大接触应力为1 954 MPa，误差为1.14%。

3 齿轮接触有限元分析

在进行齿轮接触分析时，需要考虑两种啮合状态，第1种只考虑1对轮齿啮合，第2种考虑多对齿啮合，本文目前只考虑第1种情况。根据葛宝乾[11]计算推导出的公式可知，在不考虑齿轮弯曲变形，只考虑接触变形时第1种情况接触点的载荷为

图10 有限元接触分析局部放大Fig.10 Local amplification of finite element contact analysis

(7)

式中F——接触点处的载荷

R1、R2——圆柱半径

β——主动轮转角

γ——公法线到齿轮回转中心的距离

齿轮单齿啮合最大接触应力为

(8)

ZH——节点区域系数

ZL——重合度系数

K——载荷系数

T——转矩

u——齿轮传动比

d——分度圆直径

b——齿轮厚度

设定齿轮模数m=2，齿数Z=40，压力角为20°。上面的齿轮沿着轴心逆时针方向旋转0.4°。材料的弹性模量E=2.1×105MPa，泊松比μ=0.3。

从图11～12分析结果可以看出，最大接触应力并没有出现在齿轮接触表面，而是出现在齿轮接触表面的次表层，这说明在表面切向应力的作用下，接触应力的分布发生了变化，且最大接触应力为426.3 MPa。根据齿轮单齿啮合最大接触应力公式计算可得最大接触应力为440 MPa，误差为3.1%。

图11 齿轮接触有限元应力分布Fig.11 Gear contact finite element stress distribution

图12 有限元齿轮接触分析局部放大Fig.12 Local amplification of finite element gear contact analysis

4 结论

为了能够更加精确、更加高效地得出齿轮接触分析结果，考虑了表面切向应力对接触分析的影响，考虑了四面体网格单元和六面体网格单元对分析效率的影响，并且用一个简单实例验证了这些因素的影响。从实例的分析结果可以看出：网格的类型会影响分析时的效率，四面体网格对模型分析的效率更高；表面切向应力对接触模型也有影响，表面切向应力不仅影响接触应力的大小，同时也改变了最大接触应力的分布。在以上分析得出的结论基础上，对直齿齿轮进行了有限元接触分析，得到了三维齿轮的应力分布情况和最大接触应力值，并计算出有限元分析结果与理论值之间的误差。