球与平面弹塑性接触的计算分析

官春平，金宏平

(1.广东轻工职业技术学院，广州 510300；2.湖北汽车工业学院 机械工程学院，湖北 十堰 442002)

滚动轴承作为机械设备的关键部件，其运行状态直接影响到整个机械系统的性能。据统计，约30%的旋转机械故障由滚动轴承的损伤造成。对于小型轴承，由于承受的载荷较小，通常假定轴承的内、外圈为完全刚性，而滚动体为弹性变形，然后采用Hertz接触理论研究轴承的接触应力和变形等。但对于大型轴承，如风力发电机轴承、工程机械的转盘轴承等，由于其支承结构本身受载变形较大，内、外圈滚道与滚动体在接触点会发生一定程度的塑性变形，这种情况超出了Hertz弹性接触理论的适用范围，因此必须采用弹塑性接触分析来研究轴承的接触载荷、接触应力和接触变形等。

由于轴承的结构和所承受的工作载荷均具有对称性，因此可以将其等效为单个滚动体与内、外圈的接触问题。而对于球轴承，为了便于分析，一般将其简化为球与平面的接触。根据接触力学理论，在球与平面的接触过程中，当接触变形较小时，接触体处于完全弹性接触状态，可以采用Hertz接触理论进行分析。当接触变形超过某一值后，接触体将发生屈服而出现弹塑性变形。最初的塑性变形出现在接触体内部，塑性区域位于接触中心的正上方，而接触区域仍然处于弹性变形阶段。随着外部载荷增加，接触变形进一步增大，接触体的塑性变形区域逐渐扩大，处于接触区域的材料也会发生屈服，此时，接触区域处于弹塑性变形阶段[1]。由于材料的弹塑性变形使得球与平面的接触问题复杂化，不仅存在接触的非线性和几何变形的非线性，还存在材料塑性变形的非线性，以致无法用解析模型来描述球与平面的接触过程。为此，不少学者根据接触力学和变形机理推导了多种近似解析模型，但由于目的不同，这些模型之间存在较大差别[2-6]。如果滚动体长期处于弹塑性或塑性接触状态，会严重影响轴承的寿命和旋转精度。因此，有必要对球与平面发生弹塑性接触时的参数进行研究，以便为轴承的设计或选用提供技术支持。

1 球与平面弹塑性接触理论分析

根据Tabor理论[1]，当球与平面的接触状态从完全弹性接触阶段转变到完全塑性变形阶段，即处于弹塑性变形阶段时，其最大接触应力p0将从0.6H增加到3H(H为材料的硬度)。由于该阶段接触区域既有弹性变形又有塑性变形，因此接触半径a、接触载荷F与接触变形δ之间的关系比较复杂，此时接触应力分布也不能用Hertz接触应力分布公式来描述。针对这种情况，文献[7]将接触应力进行分区假设：接触区域的中间部分存在塑性变形，其接触应力均匀分布且等于完全塑性变形时的接触应力；而接触区域的其他部分处于弹性变形阶段，其接触应力为Hertz接触应力分布，并从最大接触应力逐渐变为0。文献[8]在此基础上提出了“有限接触应力”的接触应力分布模型，如图1所示，其接触应力表达式为

图1 沿接触面的接触应力分布

(1)

式中：E1，E2分别为球与平面材料的弹性模量；ν1，ν2分别为球与平面的泊松比；Re为材料的屈服强度；R为球的半径；a为接触半径；b为接触区域中心到弹塑性接触应力分界点的径向距离。

根据接触应力连续性条件，在分界点处有

p(b)=CRe。

(2)

接触半径a远小于球半径R, 因此可得C的近似值为2.8，即在弹塑性变形阶段，接触区域中塑性变形部分的接触应力为常量。然而，根据文献[9]可知，在发生初始屈服时，p0=1.6Re。根据Tabor理论[1]，当球与平面接触处于完全塑性变形时，p0=(2.8～3)Re。因此C应该是与接触变形δ或接触半径a有关的变量，有必要对“有限应力分布”的假设进行适当修正，即

(3)

式中：k为与接触变形δ或接触半径a有关的函数。

将(3)式在整个接触区域进行积分，可以得到球与平面接触的接触载荷为

(4)

根据接触应力连续性条件，当r=b时，由(3)式可得

(5)

对(5)式进行数学变换可得

(6)

由(4)和(6)式可得

(7)

在(7)式右边，只有k和a是未知数，而k无法通过解析法求得，a无法准确测量，因而通过有限元仿真对其进行求解。

根据(7)式计算出球与平面接触的平均接触应力pm为

(8)

当球与平面的接触状态从弹性变形转变到弹塑性变形时[9]，其临界点的平均接触应力pmy为

(9)

式中：Fy和ay分别为接触变形处于临界屈服状态时的接触载荷和接触半径。

将(9)式代入(8)式，计算得到临界点处的接触应力系数k=1.466，因此k的取值范围为1.466～3。

2 有限元仿真

为提高有限元仿真的效率，减小工作量，在建模时将球与平面的接触模型简化为轴对称平面结构并采用圆柱坐标进行分析，如图2所示。在ABAQUS中采用解析刚体来模拟球。球与平面在对称轴上为轴对称约束，平面底部施加水平约束。通过对球施加z向位移来实现加载，根据球上的反作用力获得接触载荷。在仿真中，模型材料采用理想弹塑性材料。

图2 压痕模型图

3 仿真结果分析

根据Hertz接触理论[1]，球与平面弹性接触的接触半径与接触变形的关系为

(10)

式中：δy为接触处于临界屈服点时的接触变形。

球与平面的接触处于完全塑性变形阶段时，接触半径与接触变形之间的关系式为[10]

(11)

根据有限元仿真分析得到的接触半径与接触变形之间的变化关系如图3所示。从图中可以看出，当接触变形较小时，球与平面的接触处于弹性变形阶段，接触半径和接触变形之间的关系与(10)式相符。随着接触载荷的增大，接触变形增大，球与平面的接触半径也迅速增大。当接触半径a/ay=18.1时，接触面积A/Ay=2δ/δy(Ay为临界屈服时的接触面积)，此时δ/δy=163。

对图3所示的接触半径与接触变形的曲线进行拟合可得

图3 接触半径与接触变形的关系图

(12)

由此可知，当δ=δy,a/ay=0.961时，与(10)式的计算结果相比较，拟合误差为3.9%；当δ/δy=162，a/ay=17.97时，与(11)式的计算结果相比较，拟合误差为0.17%。可见，其拟合精度可以满足工程需要。

有限元仿真分析得到的最大接触应力与接触变形之间的变化关系如图4所示。从图中可以看出，当球与平面的接触处于弹性变形阶段，即δ/δy<5.4时，最大接触应力随接触变形的增大呈线性增大趋势；当δ/δy>5.4时，球与平面的接触进入到弹塑性变形阶段，接触应力增速变慢，其主要原因是接触面积快速增大。当δ/δy=163时，球与平面的接触进入到完全塑性变形阶段，此时p0/Re=2.82。

图4 最大接触应力与接触变形的关系图

对图4的曲线进行分析，并结合(3)式得到k与接触变形间的拟合函数为

(13)

平均接触应力和屈服强度之比F/(ARe)与接触变形的变化关系如图5所示。从图中可以看出，当1<δ/δy<2.55时，球与平面接触中弹性变形占主要部分，因此，可以用Hertz接触模型进行分析。当δ/δy>2.55时，球与平面接触中的塑性变形区域扩大，由于接触面积比接触载荷增大的快，因此其平均接触应力小于Hertz接触应力。

图5 平均接触应力和屈服强度之比与接触变形的关系图

根据(7)～(9)式及临界参数值，在弹塑性变形阶段有

(14)

联立(12)～(14)式，可得到弹塑性变形阶段接触载荷与接触变形间的计算关系。有限元仿真结果和Hertz接触模型计算结果对比如图6所示。从图中可以看出，当δ/δy=1时，二者误差不超过7%。

图6 仿真结果与Hertz接触模型计算结果对比图

对于半径为10 mm的GCr15钢球，其弹性模量E=210 GPa，泊松比ν=0.3，屈服强度Re=520 MPa，当其接触变形δ=0.002 7 mm时，根据(9)，(12)～(14)式，可以得到其与刚性平面发生弹塑性接触的参数，计算结果见表1。

表1 弹塑性接触参数计算结果

4 结束语

以弹性接触理论和弹塑性力学为基础，分析了球与平面发生弹塑性接触时的应力分布。采用有限元分析方法获得了弹塑性接触时的接触半径、接触应力与接触变形之间的关系。基于“有限应力分布”假设，建立了球与平面弹塑性接触的半解析模型，可获得二者间的接触半径、最大接触应力和接触载荷等参数，为球轴承的设计与选型提供参考。