CN112836445B - 一种用于求解滚动轴承的动态接触混合润滑参数的方法 - Google Patents

Abstract

本发明公开了一种用于求解滚动轴承的动态接触混合润滑参数的方法，通过建立滚动轴承模型，根据刚性位移公式和所述全局位移矩阵计算所述滚动轴承模型中每个滚动体与轴承内、外圈的接触变形量，进而计算出滚动轴承的动态载荷，进而通过动态载荷进行滚动轴承混合润滑的参数计算，以实现对轴承减摩的设计。本发明公开考虑了动态因素对滚动轴承载荷分布的影响，其计算结果更符合实际工况，同时载荷分布的计算结果直接影响混合润滑的数值计算，提高了轴承的优化性能。

Description

一种用于求解滚动轴承的动态接触混合润滑参数的方法

技术领域

本发明涉及一种滚动轴承的计算方法，具体涉及一种用于求解滚动轴承的动态接触混合润滑参数的方法。

背景技术

滚动轴承作为各类机械设备的核心部件，广泛应用于航空航天、船舶运输、铁路货运、基建等各类工程技术领域，但是目前国内轴承技术作为卡脖子技术之一，仍大量依赖进口。因此，为加快实现轴承技术的国产化，针对目前滚动轴承载荷分布及混合润滑数值计算方法相关研究的不足，提出一种用于求解滚动轴承的动态接触混合润滑参数的方法具有十分重要的意义。

在滚动轴承载荷分布研究方面，传统载荷分布计算方法认为轴承稳定运转，转速不是很高，通常忽略了动态因素的影响，例如，现有文献1[王赫男.过盈、温升和修形对城轨车辆轴箱轴承承载性能的影响[D].大连交通大学,2017.]利用有限元法对同时考虑轴与内圈过盈配合和温升影响的圆柱滚子轴承静载下载荷分布进行研究；现有文献2[付华强,杨静,陆瑞,刘桦.风电机组主轴轴承载荷分布分析及数值求解[J].风能,2019(04):120-123.]给出了风机主轴轴承静载下载荷分布的数值求解方法，并用MATLAB编译了求解程序；现有文献3[许凯,马艺宁.滚动轴承滚动体配合间隙对载荷分布的影响[J].内燃机与配件,2018(07):109-110.]采用离散方法求解静载作用下圆柱滚子轴承的载荷分布。

在轴承润滑研究方面，早期的润滑理论假设零件表面是光滑的，忽略了实际表面的粗糙峰接触，而现代工业逐渐向高速重载的方向发展，通常会选用低粘度的润滑油减小润滑油的拖拽力，但是却容易导致相互配合的工作面发生粗糙峰接触，使其工作在混合润滑状态，因此对轴承的混合润滑进行研究至关重要，目前，仅有少数学者进行相关研究，例如，现有文献4[艾路.基于粗糙表面的球轴承混合润滑研究[D].南京航空航天大学,2016.]通过渐进网格加密法(PMD)求解混合润滑问题，并讨论了不同粗糙度纹理和粗糙度参数对润滑性能的影响；现有文献5[闫晓亮.高速滚动轴承的混合润滑性能及疲劳寿命研究[D].北京理工大学,2014.]针对航空发动机主轴轴承的工况条件，研究了润滑油非牛顿特性、接触副材料、表面形貌和椭圆率等对混合润滑性能的影响。现有文献6[王骁鹏，彭文昱，田红亮，赵新泽，赵美云.基于热弹流的滤波减速器转臂轴承润滑性能分析[J].工程科学与技术，2018，50(02)：197-203.]结合点接触热弹流混合润滑理论建立轴承热弹流润滑模型，分析轴承接触副物理尺寸、载荷、卷吸速度和接触副粗糙表面等外部工况对轴承润滑特性以及表面下应力的影响。

综上所述，滚动轴承动态接触载荷分布及混合润滑的研究是高速重载轴承所关心的重点，而目前相关研究的仍然缺乏，并且轴承动态接触载荷分布直接影响混合润滑的数值计算，结合两方面进行详细的数值计算方法研究十分必要。

发明内容

本发明提出了一种用于求解滚动轴承的动态接触混合润滑参数的方法，解决了现有轴承参数优化仅考虑一方面影响因素造成优化结果不准确的问题。

本发明采用的技术手段如下：

一种用于求解滚动轴承的动态接触混合润滑参数的方法，包括以下步骤：

步骤1、建立滚动轴承模型，给定所述滚动轴承模型的全局位移矩阵的初始值，通过刚性位移公式和所述全局位移矩阵计算所述滚动轴承模型中每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o；

步骤2、通过所述每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o，计算每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o；

步骤3、通过所述每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o计算局部载荷和；

步骤4、计算所述局部载荷和与实际工况外载荷差值，判断结果是否收敛，若是，输出所述每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o，并执行步骤5；若否，更新所述全局位移矩阵，并重新执行步骤1；

步骤5、建立滚动轴承混合润滑统计学模型，并输入初始参数，所述初始参数包括每个滚动体与轴承内圈的接触载荷Q_i、每个滚动体与轴承外圈的接触载荷Q_o、速度、材料特征和表面特征参数，设置压力分布和最小膜厚的初值；

步骤6、通过所述滚动轴承混合润滑统计学模型计算润滑油粘度、密度、表面弹性变形和膜厚分布；

步骤7、根据计算的所述润滑油粘度、密度、表面弹性变形和膜厚分布并通过求解离散化的平均Reynolds方程和KE粗糙峰接触公式计算压力分布；

步骤8、判断压力是否收敛，若收敛则进行步骤9，若不收敛则修正膜厚后返回步骤6；

步骤9、通过求解离散化的能量方程计算温度分布；

步骤10、判断温度是否收敛，若收敛则输出膜厚、温度、粗糙峰摩擦阻力、流体摩擦阻力以及摩擦阻力矩以实现对轴承性能评价，若不收敛则修正温度后返回步骤6。

与现有技术比较，本发明公开的用于求解滚动轴承的动态接触混合润滑参数的方法具有以下优点，1、本公开考虑了动态因素对滚动轴承载荷分布的影响，其计算结果更符合实际工况，同时载荷分布的计算结果直接影响混合润滑的数值计算；2、本公开以动态接触载荷分布作为初始条件之一，提出了一种混合润滑的数值计算方法，该方法充分考虑载荷、转速、游隙及零件表面粗糙度对轴承润滑的影响；3、本公开以MALTLB软件作为数值计算工具，大大提升计算效率。

附图说明

图1是本发明的滚动轴承动态接触载荷分布计算的流程图；

图2是本发明的滚动轴承动态接触载荷分布计算的整体思路流程图；

图3是本发明的滚动轴承动态接触载荷分布计算的收敛判断流程图；

图4是本发明的滚动轴承混合润滑计算的流程图；

图5是实施例1中6016深沟球轴承模型图；

图6是实施例1中6016深沟球轴承各编号滚动体位置示意图；

图7是实施例1中6016深沟球轴承中滚子与内圈接触的等效模型图；

图8是实施例2中的等效模型及几何尺寸图；

图9是实施例2中的模型工作条件的MATLAB设定程序。

具体实施方式

如图1、图2、图3和图4所示，一种用于求解滚动轴承的动态接触混合润滑参数的方法，包括以下步骤：

步骤1、建立滚动轴承模型，给定所述滚动轴承模型的全局位移矩阵的初始值，通过刚性位移公式和所述全局位移矩阵计算所述滚动轴承模型中每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o；

步骤2、通过所述每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o，计算每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o；

具体地，将每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o代入Hertz公式δ＝KQ^2/3，可以计算每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o；

步骤3、通过所述每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o计算局部载荷和；

具体地，该步骤中通过全局受力平衡方程计算局部载荷和，全局受力平衡方程为动力学微分方程组，具体如下：

步骤4、计算所述局部载荷和与实际工况外载荷差值，判断结果是否收敛，若是，输出每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o，并执行步骤5；若否，更新所述全局位移矩阵，并重新执行步骤1；

具体地，如图3所示，对轴承整体施加位移d0后，开始计算轴承各个滚动体所受的载荷及变形，这一步是利用变形确定的载荷，计算载荷过程中需要计算刚度K，然后把计算的各个滚动体载荷求和得到fn，若各个滚动体载荷求和得到fn与实际外载f相对误差小于设定值，本实施例中相对误差设为10e-6，则认为收敛，输出各个滚动体载荷及变形，若不收敛，调整初始位移继续迭代，直到计算收敛。

步骤5、建立滚动轴承混合润滑统计学模型，并输入初始参数，所述初始参数包括每个滚动体与轴承内圈的接触载荷Q_i、每个滚动体与轴承外圈的接触载荷Q_o、速度、材料特征和表面特征参数，设置压力分布和最小膜厚的初值；

具体的，滚动轴承混合润滑统计学模型中含有大量控制方程，其控制方程包括：平均Reynolds(雷诺)方程、能量方程、粗糙峰接触公式、弹性变形方程、膜厚方程、载荷平衡方程及润滑油粘压-粘温及密压方程；

(1)其中，压力大小通过平均Reynolds方程求解，当卷吸速度方向(即短轴)与x轴方向相同时，方程可简化为如下形式：

式中，u_s＝|u₁-u₂|为两表面的相对运动速度，u₁、u₂分别为两表面的速度。

(2)求解膜厚时考虑表面粗糙度的影响，具体膜厚方程为：

其中，v(x，y.t)为弹性变形：

式中，h₀为刚体中心膜厚，E是上下接触面的等效弹性模量，R_x、R_y是两接触在x轴和y轴两个方向上综合曲率半径，p为流体润滑压力，δ(x，y，t)为确定型粗糙函数。

(3)载荷平衡方程为：

∫∫p(x，y，t)＝w

(4)粘压-粘温方程采用Roelands粘压-粘温关系式，具体为：

式中：Z为粘压系数，S为粘温系数，T₀为初始温度，一般取303K。

(5)密压方程为：

(6)能量方程为：

上下界面温度条件：

式中，k为润滑油热传导系数，k₁，k₂分别为上下界面热传导系数，ρ₁，ρ₂分别为上下表面的宽度，c₁，c₂分别为上下表面的比热系数，u₁，u₂分别为上下表面的线速度，T₀为初始温度。

(7)粗糙峰经历弹性变形、弹塑性变形到完全塑性变形的全过程，粗糙峰的接触压力可表示为：

式中：

带星号的变量表示无量纲化，被除因子为σ；nβσ＝0.05，β为粗糙峰的半径，n为粗糙峰的密度；HV为较软材料的维氏硬度；d₁为开始产生屈服点的位移；d₂为产生完全塑性流动的位移，且d₂＝54d₁；σ_s为表面粗糙度峰值点的均方根；y_s为表面高度中线与表面峰值点中线的距离。采用有限差分法将混合润滑统计学模型的上述全部控制方程离散为一个线性方程组，通过数值迭代求解压力分布、膜厚分布及温度分布，并利用MALTLB软件进行数值计算；

步骤6、通过所述滚动轴承混合润滑统计学模型计算润滑油粘度、密度、表面弹性变形和膜厚分布；

步骤7、根据计算的所述润滑油粘度、密度、表面弹性变形和膜厚分布并通过求解离散化的平均Reynolds方程和KE粗糙峰接触公式计算压力分布；

步骤8、判断压力是否收敛，若收敛则进行步骤9，若不收敛则修正膜厚后返回步骤6；

步骤9、通过求解离散化的能量方程计算温度分布；

步骤10、判断温度是否收敛，若收敛则输出计算结果以实现对轴承的性能评价，若不收敛则修正温度后返回步骤6。计算结果包括最小膜厚、中心膜厚、膜厚比、温度以及摩擦阻力，通过这些参数可以实现对轴承的性能评价，进而可以进一步用于对轴承进行优化。

以下结合附图和实施例对本发明做进一步说明，分别选取点接触的6016深沟球轴承及线接触的NU210圆柱滚子轴承作为实施例。

实施例1：如附图5所示，6016深沟球轴承共有14个钢球，根据计算边界条件显示，此轴承既承受轴向载荷(Fa＝1019.9N)，同时也承受径向载荷(Fr＝3403.45N)，钢球的接触数为14个，根据6016深沟球轴承在实际工况下的边界条件及轴承三维模型，分别计算轴承在额定转速(5433rpm/min)以及最高转速(7960rpm/min)情况下，轴承载荷分布、膜厚、摩擦阻力等参数随速度的变化。本实施例还考虑了轴承游隙(80um、95um、105um)带来的影响，游隙的不同导致滚动体接触角发生了变化，最终引起滚动体载荷的略微差异。

对实施例1进行以下步骤，直至结果收敛。

步骤(1)：给予滚动轴承模型一个全局位移初始矩阵d₀，通过刚体位移公式计算各滚动体变形量δ_i与δ_o；

步骤(2)：将变形量代入Hertz公式δ＝KQ^2/3，计算得到每一个滚动体承受的载荷Q_i与Q_o；

步骤(3)：利用全局受力平衡方程，判断结果是否收敛(|δ|＜10^-6)，如不收敛则进行迭代循环计算直至收敛。

实施例1中各编号滚动体位置示意图如附图6所示，不同转速及游隙下实施例1轴承载荷分布结果如表1。

表1不同转速下各滚动体载荷(N)

混合润滑计算流程如附图4所示，具体分为以下步骤：

步骤(1)：首先输入初始参数，载荷、速度、材料特征和表面特征参数等，合理设置压力分布和最小膜厚的初值；

步骤(2)：计算润滑油粘度、密度、表面弹性变形和膜厚分布等；

步骤(3)：通过求解离散化的平均Reynolds方程和KE粗糙峰接触公式更新压力分布；

步骤(4)：判断压力是否收敛，若收敛则进行下一步，若不收敛则重新进行压力、膜厚等迭代计算随后通过求解离散化的能量方程计算温度温度分布最后判断温度是否收敛，若收敛则输出计算结果，结束程序，若不收敛则重新返回进行压力和膜厚的迭代计算。

实施例1中选取最小膜厚、中心膜厚和摩擦阻力矩为混合润滑计算结果。

(1)膜厚

最小膜厚采用Matlab语句Hmin＝min(min(H))*aB^2/Rx，计算出单个滚动体模型收敛后，膜厚H的最小值。其中，H为量纲一化的膜厚分布，aB^2/Rx将膜厚转变为实际值。

中心膜厚采用Matlab语Hcen＝H((N+1)/2，(N+1)/2)*aB^2/Rx计算出单个模型收敛之后，中心点膜厚Hcen的值。其中，H((N+1)/2，(N+1)/2)为XOY坐标系下中心点的膜厚无量纲值，aB^2/Rx将膜厚转变为实际值。

(2)摩擦阻力矩

摩擦阻力矩方面，首先利用动平衡下轴承载荷分布计算单个滚动体所承受的载荷，接着利用模型计算得出滚动体与内滚道之间的摩擦力，包括流体的摩擦阻力以及接触粗糙峰的摩擦阻力两部分，接着计算出摩擦阻力矩，最终将承载滚动体的力矩进行叠加即可得到所有承载滚动体和内圈之间的摩擦所产生的总体摩擦扭矩。

流体的摩擦阻力(内圈部分)为：

接触粗糙峰的摩擦阻力(内圈部分)为：

F_contact＝∫∫μ_cp_cdxdy

单个滚动体与内滚道的摩擦阻力为：

F_total＝F_fluid+F_contact

单个滚动体与内滚道的摩擦阻力矩为：

T＝F_total*R

所有滚动体与内滚道产生的总摩擦阻力矩为：

实施例1轴承摩擦阻力矩由于对称原因，计算公式如下

M＝∑T_i＝R[F_NO.1+F_NO.8+2(F_NO.2+F_NO.3+F_NO.4+F_NO.5+F_NO.6+F_NO.7)]

(3)温度

温度方面，已知润滑油环境温度为80～100℃，因此设定本计算模型初始温度为90℃。利用MATLAB语句，导出收敛后各滚动体最高温度，和中心点温度。

实施例1游隙80um额定转速下(5433rpm/min)，轴承混合润滑各参数计算结果如表2所示：

表2游隙80um轴承混合润滑各参数计算结果

实施例2：牵引电机中的NU210圆柱滚子轴承，滚子尺寸为11mm*12mm，内圈直径为50mm，外圈直径为90mm，在三种不同工况(FR＝2.645KN、4.978KN、6.541KN)进行额定转速及最高转速下的计算分析。

实施例2为有限长线接触模型，由于接触区的宽度远小于接触点的曲率半径，因而可以对接触表面做适当的几何简化。本模型将圆柱滚子和内圈等效为一个当量圆柱与一个平面的接触，如附图7和图8所示为等效模型及模型几何尺寸。

实施例2动态接触载荷分布的计算方法与实施例1相同，利用MATLAB软件建立圆柱滚子轴承载荷分布的程序，其中模型工作条件的设定如附图9所示，计算出的NU210轴承在三种不同工况下的载荷分布情况如下表所示。

表3三种工况下的载荷分布计算结果

实施例2混合润滑计算步骤与实施例1基本相同，同样选取最小膜厚、中心膜厚和摩擦阻力矩为计算结果。

(1)膜厚

最小膜厚采用Matlab语句Hmin＝min(min(H))*B^2/Rx，计算出单个滚子模型收敛后，膜厚H的最小值。其中，H为量纲一化的膜厚分布，B^2/Rx将膜厚转变为实际值。

中心膜厚采用Matlab语Hcen＝H((XN+1)/2，(YNN+1)/2)*B^2/Rx计算出单个模型收敛之后，中心点膜厚Hcen的值。其中，H((XN+1)/2，(YNN+1)/2)为XOY坐标系下中心点的膜厚无量纲值，B^2/Rx将膜厚转变为实际值。

(2)摩擦阻力矩

摩擦阻力矩方面，首先利用上述公式计算单个滚动体所承受的载荷，接着利用模型计算得出滚动体与内滚道之间的摩擦力包括流体的摩擦阻力以及接触粗糙峰的摩擦阻力两部分，接着计算出摩擦阻力矩，最终将承载滚子的力矩进行叠加即可得到所有承载滚子和内圈之间的摩擦所产生的总体摩擦扭矩。

流体的摩擦阻力(内圈部分)为：

接触粗糙峰的摩擦阻力(内圈部分)为：

F_contact＝∫∫μ_cp_cdxdy

单个滚动体与内滚道的摩擦阻力为：

F_total＝F_fluid+F_contact

单个滚动体与内滚道的摩擦阻力矩为：

T＝F_total*R

所有滚动体与内滚道产生的总摩擦阻力矩为：

实施例2工况一(FR＝2.645KN)时，轴承混合润滑各参数计算结果如表4所示：

表4工况一(FR＝2.645KN)轴承混合润滑各参数计算结果

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种用于求解滚动轴承的动态接触混合润滑参数的方法，其特征在于：包括以下步骤：

步骤1、建立滚动轴承模型，给定所述滚动轴承模型的全局位移矩阵的初始值，通过刚性位移公式和所述全局位移矩阵计算所述滚动轴承模型中每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o；

步骤2、通过所述每个滚动体与轴承内圈的接触变形量δ_i和每个滚动体与轴承外圈的接触变形量δ_o，计算每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o；

步骤3、通过所述每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o计算局部载荷和；

步骤4、计算所述局部载荷和与实际工况外载荷差值，判断结果是否收敛，若是，输出所述每个滚动体与轴承内圈的接触载荷Q_i和每个滚动体与轴承外圈的接触载荷Q_o，并执行步骤5；若否，更新所述全局位移矩阵，并重新执行步骤1；

步骤5、建立滚动轴承混合润滑统计学模型，并输入初始参数，所述初始参数包括每个滚动体与轴承内圈的接触载荷Q_i、每个滚动体与轴承外圈的接触载荷Q_o、速度、材料特征和表面特征参数，设置压力分布和最小膜厚的初值；

步骤6、通过所述滚动轴承混合润滑统计学模型计算润滑油粘度、密度、表面弹性变形和膜厚分布；

步骤7、根据计算获得的所述润滑油粘度、密度、表面弹性变形和膜厚分布并通过求解离散化的平均Reynolds方程和KE粗糙峰接触公式计算压力分布；

步骤8、判断压力是否收敛，若收敛则进行步骤9，若不收敛则对膜厚进行修正后返回步骤6；

步骤9、通过求解离散化的能量方程计算温度分布；

步骤10、判断温度是否收敛，若收敛则输出膜厚、温度、粗糙峰摩擦阻力、流体摩擦阻力以及摩擦阻力矩以实现对轴承性能评价，若不收敛则修正温度后返回步骤6。

Images