APP下载

圆柱滚子轴承数值求解方法的研究与分析

2023-02-02李猛帅郑甲红

科技与创新 2023年2期
关键词:柔度滚子单元格

李猛帅,郑甲红

(陕西科技大学,陕西 西安 710021)

1 确定相关参数

圆柱滚子轴承的接触问题属于有限长线接触问题,此类问题是在基于点接触问题的基础上展开讨论的,因此点接触问题的一些结论对线接触是有用的。包括以下参数:主曲率半径R、变形前后滚子表面对应点与起始接触点之间的距离Z、接触应力与外载荷之间的平衡方程。

主曲率半径是根据圆柱滚子轴承的实际滚子状态来确定的,变形前后滚子表面对应点与起始接触点之间的距离与点接触不同,它在计算程序中作为初始参数输入,而且需要将其在计算程序中转化成向量的形式,它的确定方式如下。

采用文献[1]中的式1.21—式1.22b得出:

式(2)(3)中:R11、R12、R21、R22分别为2个接触物体在一点上的主曲率半径。

对于圆柱滚子轴承而言,接触角α=0°,R12、R21、R22都为无穷大,代入式(2)与式(3)并将2式相加可得:

这样就确定了Z的表达式,在MATLAB软件中输入时,需要将其以向量的形式输入,而且它的个数也要和2种数值求解方法中的网格数一一对应,因为这样才能满足矩阵之间的基本运算,从而进一步求解基本平衡方程。

线接触问题的基本平衡方程与点接触问题的基本方程在原理上是类似的,可理解为是将点接触问题的基本方程转化成求和的形式。

2 非Hertz接触问题的求解方法

有限长线接触问题中短圆柱滚子的接触状态一般是在非理想状态下的,实际接触状态比较复杂。一般在求解时,是将所有可能接触区域划分成多个矩形小单元格,而且要注意所求解区域要稍微大于实际接触区域,这样也符合实际工程要求。

划分矩形单元如图1所示,将可能接触的区域划分成多个矩形小单元,其中任意一个单元j的半长和半宽分别为aj和bj,假设单元上的接触应力为pj[1]。

图1 划分矩形单元

则单元j上的接触应力在任意一个单元i的中心上产生的位移为:

式(6)中:Dij为柔度系数。

柔度系数Dij通过MATLAB软件进行计算,而且它的下标在MATLAB中代表着任意一个矩阵,也就是单元矩阵划分多少个,那么对应的柔度系数矩阵在MATLAB中就有多少个。那么线接触的基本方程就可以表示为:

式(7)(8)中:m、n为单元格数目,它们与柔度系数矩阵下标表示的意义是一样的,在MATLAB中输入它们时,将二者设定为相等数目;E′为弹性模量;δ为弹性趋近量。

需要注意的是,在计算时要满足pj≥0,因为当pj<0时,代表此单元格已经超出实际接触区域,没有受到外载荷的作用,在下一次程序迭代时就不再进行计算。

整个关键参数在MATLAB程序中运算时,按以下流程进行:首先输入外载荷Q以及物理几何参数,计算向量Zi,并且将柔度系数Dij计算出来后转化成矩阵形式;程序中δ虽然作为未知量,但是它是通过基本方程和误差修正最终确定的,因此需要赋初值;其次通过解线接触的基本方程求出接触应力pj,并判断pj是否小于0,若pj>0则进行下一步,如果pj<0,则将pj的值赋为0,循环到解式(8);在解方程时,pj=0表示在接下来的计算时,将pj<0的单元格跳过,直接进行下一单元的计算,以此类推。最后一步的判断条件为取ε为0.001,若满足条件,则结束程序;若不满足,循环跳到求解式(8),重复进行循环,直至将所有的接触区域解算完毕。

3 Hertz接触问题的求解方法

利用Hertz理论同样也可以求解线接触问题,式(7)和式(8)是接触问题积分方程数值解的最简单形式之一,它已被许多学者以不同的形式成功地应用于分析Hertz问题和非Hertz问题[2-5]。然而基于对线接触问题的理解,还可以对圆柱滚子轴承有限长线接触问题的数值解进行进一步的简化,这也就是对圆柱滚子有限长线接触问题进行处理的一维处理方法。

它的求解思想也是将所有可能的接触区域进行单元格划分,不同的是在单元j内假定接触应力沿滚子素线方向(y轴)为均匀分布,沿横向(x轴)按Hertz分布[1],条形单元划分如图2所示。

图2 条形单元划分

图2中p0j为单元中心处的最大接触应力,单元j上的应力在其他任意单元格i的中心处产生的位移与非Hertz接触一样。

式(9)中:hj为单元格宽度的1/2;aj为半长。

此处柔度系数与Hertz不同,主要是在x轴方向上,将应力分布看成是按Hertz分布的,那么就可以先对y′积分,然后对x′积分,就可以求解出柔度系数,再根据单元格划分,将Dij在MATLAB软件中输入,并且同样转化为矩阵的形式。

Hertz线接触的基本方程为:

它也是将点接触的方程转化为求和的形式,而且与非Hertz方程比较相似。

此套计算程序比非Hertz计算稍微复杂,主要流程为:首先输入外载荷Q、向量Zi、网格划分个数n及几何材料参数,并且对δ赋初值,求解向量Si=[δ-Zi],这一步是为了简化式(11)的计算,并且对aj赋初值;其次计算柔度矩阵,求解式(11),并判断解出的p0j是否大于0,大于0则进行下一步,小于0时回到计算柔度矩阵这一步重新计算;如果满足条件求解式(10),若不满足则同样跳回到计算柔度系数矩阵这一步,重新循环,目的是为了更精准地自动跳过非接触区域的计算;最后,与非Hertz问题一样,判断根据需要调整δ的值,从程的值,从程序第一步重新计算。

4 实例计算

某二维转台的圆柱滚子轴承,已知滚子曲率半径R11=10 mm,R11=R21=R22=∞,外载荷Q=2 000 N,单元格划分个数n=18,弹性模量E=2.06×105N/mm2,滚子素线长L=30 mm。计算该圆柱滚子在外载荷作用下的接触应力及弹性趋近量。

根据编制好的2种计算程序,将所有初始条件输入到MATLAB软件中,通过对2种计算方法的仔细分析,把对结果有直接影响的参数一一列举,并整合出曲线图。网格划分个数与单元格中心最大接触应力之间的关系如图3所示。

图3 网格划分与最大接触应力之间的关系

由图3可知,利用2种计算方法进行单元格划分时,不同单元格中心点处的最大接触应力是不同的,而且它们都随着网格数的增加整体呈逐渐递增的趋势。两者的数值比较接近,由于在非Hertz计算程序中,它的网格划分数不能超过27个,当网格个数超过27个时,就会导致柔度系数矩阵的行列式为0,程序就无法进行计算,因此具有一定的局限性,而Hertz计算程序中网格划分数可以取到更多。

圆柱滚子的主曲率半径与弹性趋近量之间关系的曲线图如图4所示。由图4可知,随着圆柱滚子的主曲率半径逐渐增大时,利用2种计算方法得出的弹性趋近量都呈逐渐递减的趋势,但是在非Hertz计算方法中曲线下降比较明显,而Hertz计算中曲线下滑比较缓慢,弹性趋近量的随主曲率半径增大波动很小,2种计算方法虽然有一定的误差,但是误差较小。一般采用Hertz计算方法进行求解。

图4 主曲率半径与弹性趋近量之间的关系

网格划分个数与弹性趋近量之间关系的曲线图如图5所示。由图5可知,网格划个分数不同时,利用2种计算方法中得出的弹性趋近量也是有差异的,随着矩形单元格个数增加,弹性趋近量都呈递增趋势,而且2种计算方法的结果比较接近,这里也需要注意在非Hertz计算方法中单元格的个数不能超过27个,否则会影响计算结果。

图5 网格划分与弹性趋近量之间的关系

主曲率半径与单元格中心最大接触应力之间的关系如图6所示。由图6可知,主曲率半径改变时,2种计算方法中的单元格中心处的最大接触应力也随着改变,随着主曲率半径逐渐增大而呈明显上升趋势,当主曲率半径为70~80 mm之间时,最大接触应力有下降趋势,在80 mm之后又逐渐增大。

图6 主曲率半径与单元格中心最大接触应力之间的关系

5 结论

弹性趋近量和接触应力是后期求解轴承刚度的基础,在2种计算方法中,对结果有影响的主要参数包括主曲率半径和网格划分个数,随着网格划分个数逐渐增多,对应的单元格中心处的最大接触应力逐渐增大,虽然有个别突然减小,但整体上都呈上升趋势;对应的弹性趋近量也随着单元格个数增多而逐渐增大。当主曲率半径逐渐增大时,两者的弹性趋近量都逐渐减小,其中非Hertz计算方法减小比较明显,而Hertz计算方法中减小幅度较小,2种曲线虽然看上去差别较大,但数值相差较小。2种计算方法虽然存在一定的误差,但误差较小,可以满足实际要求。主曲率半径对结果的影响在后期轴承的选择上有一定的借鉴意义,根据需要对网格划分个数进行合理的划分和求解。

猜你喜欢

柔度滚子单元格
自重荷载下非均匀支撑板式无砟轨道静态响应
特大型调心滚子硬车削工艺试验探究
圆锥滚子轴承半凸滚子的优化研究
流水账分类统计巧实现
仿真模拟在多联推力滚子轴承研发中的应用
圆柱滚子轴承失效分析
基于ANSYS的新型椭圆铰链疲劳仿真分析
玩转方格
玩转方格
新型直圆导角复合型多轴柔性铰链的柔度计算及其性能分析