滚动轴承运动-接触精确建模与二次开发

马帅军,闫柯,张晓红,王明凯,朱永生,洪军

(西安交通大学现代设计及转子轴承系统教育部重点实验室,710049,西安)

滚动轴承作为旋转机械转子系统的核心支撑零部件,其动力学性能直接影响装备整机的运行状态。轴承力学的发展经历了静力学、拟静力学、拟动力学和动力学四个阶段,其中动力学阶段考虑因素最多,计算精度最高,最为完善[1]。然而,基于运动微分方程对轴承动力学特性进行描述时,需要对轴承各组件进行复杂的空间力系分析,且存在方程参数多、初值依赖性强、非线性弱、收敛性不强等问题,例如广为所知的动力学软件ADORE[2],在仿真过程中需要设置配合公差、摩擦拖动系数、疲劳系数等上百个参数,操作专业性极强,给工程应用带来一定不便。随着多体动力学软件的不断成熟,如ADAMS、Simpack、RecurDyn等,通过将轴承各个零件视为纯刚体,忽略大变形,细化小变形,大幅度提升了动力学分析软件的计算效率[3],成为当前滚动轴承工程应用分析中的首选方案。

国内学者基于ADAMS等软件,围绕滚动轴承动态特性开展了大量研究工作。关猛等利用Pro-Engineer软件建立双列圆锥滚子轴承的三维模型,导入到ADAMS软件进行相关参数设置,进而开展了轴承的动力学特性仿真分析[4];姚廷强等直接利用ADAMS软件中的CMD参数化命令建立了四点角接触球轴承三维模型,并利用宏命令进行了相关参数设置,分析了不同工况下轴承的动态性能[5];张军飞等对上述仿真步骤进行了集成,编写了简单的轴承仿真界面,包括参数化建模和参数设置功能,提升了轴承动力学分析的便利性[6]。考虑到滚动轴承组件间的点线高副接触特征,其接触变形直接影响着套圈与滚子之间的接触载荷与油膜压力分布特点,进而影响轴承磨损、发热、疲劳等一系列服役性能预测的准确性[7],因此,在利用ADAMS软件开展轴承动力学仿真时,解决轴承套圈与滚动体间的接触问题显得尤为重要。ADAMS软件在处理接触问题时,有两种库可供选择:Default_Library和Parasoilds,两者之间的主要区别如图1所示。其中Dafault_Library对体的识别精度较低,适合分析宏观接触,但计算效率高;Parasoilds对体的识别精度较高,适合分析微观接触,但计算效率低,很难在工程中应用[8]。

(a)Default_Library库

(b)Parasoilds库

针对上述问题,李琦利用Fortran语言编写了Gfosub作用力子程序,对ADAMS进行了二次开发,开展了轴承动力学性能研究,但仍然没有完全脱离软件自身的接触设置,仿真分析精度较低[9];邓四二等利用Fortran语言编写了Gfosub作用力子程序,对ADAMS进行了二次开发,没有依赖于软件自身接触,计算精度较高[10],但是建立的动力学模型是基于局部坐标系的,在处理形状较为规则的保持架与滚动体时效果较好[11],而对于形状复杂的异型保持架、修型滚子等稍显不足。

针对ADAMS软件自身接触求解器的问题,本文基于轴承全局坐标系,利用C语言编写Gfosub作用力子程序,以角接触球轴承为例,建立了轴承动力学模型,编写的子程序与软件自身的数据传递原理如图2所示,针对Gfosub子程序及ADAMS求解器的特点,选取GSTIFF I3求解方法来求解建立的动力学模型。通过将所建立模型的仿真结果与ADAMS自身接触模型、经典的轴承力学模型结果对比,证明了本模型的接触仿真计算精度较高。最后搭建保持架转速测量实验台,验证本模型的运动仿真分析精度。

1 动力学模型的建立

角接触球轴承动力学模型主要由各个零件之间的相对位置关系以及滚动体、保持架、内外套圈之间的作用力构成。对模型作如下假设。

(1)轴承各个零件的质心与其几何中心重合。

(2)轴承各个零件均为刚体,不考虑轴承零件的大变形。

图2 ADAMS二次开发原理示意图Fig.2 A schematic diagram of ADAMS secondary development

(3)外圈固定,限制6个自由度;内圈由于需要添加驱动转速,故具有5个自由度;滚动体具有6个自由度;保持架具有3个自由度。

1.1 轴承中各个参考坐标系的建立

为了能够准确地确定轴承各个零件之间的几何位置关系,建立角接触球轴承的整体坐标系与各个局部坐标系。各零件坐标系如图3所示。

图3 角接触球轴承的各零件坐标系Fig.3 Coordinate system for each part of angular contact ball bearing

(1)Oxyz为轴承整体的惯性坐标系,O点与外圈的质心相重合(假设外圈固定,内圈旋转),x轴为轴承的轴向,y轴垂直于x轴指向1号滚动体,z轴方向可根据右手定则确定。

(2)Oixiyizi为内圈的随体坐标系,Oi点与内圈的质心重合。

(3)Oajxajyajzaj为滚动体的随体坐标系,Oaj点与滚动体j的质心重合。

(4)Objxbjybjzbj为滚动体的方位坐标系,Obj点与滚动体j的质心重合,xbj轴时刻指向轴承的轴向,ybj轴垂直于xbj径向外指,zbj轴由右手定则确定其方向。

(5)Ocxcyczc为保持架的中心坐标系,Oc点与保持架整体的形心重合。

(6)Ocjxcjycjzcj为保持架的兜孔坐标系,Ocj点与兜孔中心重合。

1.2 轴承各个零件之间的相互作用力

(1)

式中:Toi为内圈连体坐标系到轴承整体坐标系的旋转变换矩阵;r2为滚动体中心与内圈中心的距离;ϑj为滚动体j相对于轴承内圈坐标系下的方位角。

(2)

式中:ri为内圈的沟曲率半径;Dw为滚动体直径。

图4 各个零件之间的相互位置关系Fig.4 Position relationship between parts

虽然滚动体与滚道之间被弹流油膜分离,但对其接触参数影响不大,仍可用Hertz理论进行接触力的计算。第j个滚动体与内滚道之间的法向接触力为

(3)

式中:Kij为第j个滚动体与内滚道间的Hertz接触刚度,具体计算详见文献[12]。

同理可得第j个滚动体与外滚道之间的法向接触力

(4)

1.2.2 滚动体与套圈之间的摩擦力 在滚动体与内圈滚道的椭圆接触面上,由于接触点处两者的滑动速度不同,故存在相对运动趋势或者相对运动,导致滚动体与内外滚道表面产生油膜拖动力。第j个滚动体与内圈滚道表面的拖动力Tix(y)j可由滚动体与滚道之间的接触应力和拖动系数的乘积进行积分得到,即为

(5)

同理得第j个滚动体与外滚道之间的拖动力为

(6)

(7)

式中:Dcm为保持架兜孔中心所在的圆的直径。

图5 保持架与滚动体之间的几何关系Fig.5 Geometric relationship between cage and roller

(8)

式中:Cp为保持架与滚动体之间的兜孔间隙。

保持架的兜孔和滚动体接触的模型可视为Hertz接触,接触刚度Kcj可参照套圈与滚动体的接触刚度计算方法,则两者之间的接触力为

(9)

1.2.4 保持架与引导套圈之间的相互作用力 保持架与引导套圈之间的相互作用是由于润滑剂的流体动压效应所产生的,而套圈引导面和保持架端面的作用力比较小,故可采用短滑动轴承模型[15]。保持架与引导套圈的相关几何关系如图6所示。

图6 保持架与引导套圈之间受力示意Fig.6 Force diagram between cage and guide ring

在图6中,平面坐标系Sc=(Oc,yc,zc)固定在保持架上,yc轴指向由于流体动压效应产生的最小油膜厚度h0所在点。以外引导的保持架为例,某一时刻保持架所受的力与力矩为

(10)

式中:η0为一般环境下润滑油的动力黏度;u1为润滑油的拖动速度;R1、L分别为保持架定心表面半径和宽度;C1为保持架的引导间隙;ε为保持架质心与中心的偏移量。

为了较为方便地建立保持架的动力学微分方程,将保持架所受的力转换到轴承的整体坐标系中,其作用力可表示为

(11)

式中:Tc为保持架坐标系到整体坐标系的变换矩阵。

1.2.5 滚动体的润滑黏性阻力 在轴承运行过程中,轴承腔内会存在大量的润滑剂以及空气组成的混合物,容易对滚动体产生相关的黏性阻力,对轴承性能产生影响。但是由于内部几何形状、运行环境都比较复杂,建立准确的黏性阻力计算模型比较困难。在此将滚动体受到的黏性阻力等同于圆柱体在流体中移动所受的阻力[16-17],故可按下式计算:

(12)

式中:wmj为滚动体公转速度;Cd为圆柱体的黏性阻力系数;ρ0为润滑油密度;dm为轴承的节圆直径;S为圆柱体的平移方向面积。

1.3 微分方程的建立

(1)对保持架进行受力分析,如图7所示,即可得到保持架的运动微分方程

(13)

式中:mc为保持架质量;Ic为保持架转动惯量;yc、zc为保持架位移;wcx为保持架公转角速度。

图7 某时刻保持架的受力情况Fig.7 The force of the cage at a certain time

(2)对内圈进行受力分析,如图8所示,即可得到内圈的运动微分方程

(14)

图8 某时刻内圈的受力情况Fig.8 The force on the inner ring at a certain moment

式中:mi为内圈的质量;FX、FY、FZ、MY、MZ为施加在内圈上的载荷;Iiy、Iiz为内圈在各个方向的转动惯量。

2 动力学模型在动力学软件中的实现

2.1 作用力子程序的开发

ADAMS软件不仅提供了易于操作的用户界面以及强大的动力学求解功能,还留有特定的接口,方便用户针对特殊需求基于软件进行二次开发。然而轴承各个零件之间的受力较为复杂,需要以子程序的形式来定义零件之间的相互作用力。在开发轴承作用力子程序时,采用ADAMS软件提供的用户作用力子程序的模板进行编写,这个子程序需要包含一个头文件slv_c_utils.h,用于预定义在子程序中用于子程序和ADAMS软件自身的数据传输函数[18-19]。在进行子程序编写时,仿真过程中调用的作用力子程序中的函数名称、形参个数和类型必须按照ADAMS预先规定的格式进行编写。在子程序中通过SYSARY函数或SYSFNC函数提取求解器中的状态变量,并且利用RESULT数组函数将计算的结果返回到求解器之中。

2.2 求解器的选择

基于Gfosub编写完成作用力子程序之后,需要在ADAMS软件建立微分方程,并完成求解工作。而在ADAMS求解器对微分方程进行求解时,主要包含以下3个步骤:

(1)利用显式求解方法对微分方程进行初步求解,得到微分方程解的预估值;

(2)利用隐式求解方法再次对微分方程进行求解,得到微分方程解的校正值;

(3)将预估值与校正值之间的误差与求解器设定误差进行对比,若满足设定的误差要求,则作为微分方程的解,若不满足,则缩小步长并返回步骤1,继续求解。

在求解过程中,ADAMS一共提供5种类型的预估-校正方法:GSTIFF、WSTIFF、HASTIFF、HHT、Newmark,各种方法的具体特点如表1所示。其中GSTIFF和WSTIFF的I3积分格式求解器的特点是将误差收敛在位移上,GSTIFF、WSTIFF的I2积分格式和HASTIFF的SI1/SI2求解器的特点是将误差收敛在速度上,HHT和Newmark求解器是将误差收敛在加速度上[20]。基于轴承的相关基础理论,第一种仿真计算精度最高,第三种仿真计算精度最差,在此选择GSTIFF和WSTIFF的I3积分格式求解器。

表1 ADAMS中各种求解方法的特点

2.3 Gfosub子程序与自身接触设置对比

为了验证本文所提出的建模方法的优越性,与常规直接应用ADAMS软件自身的接触算法(Dafault_Library库)进行求解的模型进行对比分析。以7008C轴承为例,利用两种方法建立相应的仿真模型,设置相同的仿真工况:ni=4 000 r/min,Fa=450 N,Fr=500 N。选取模型中1号滚动体与内、外圈的接触力计算结果进行对比,结果如图9和图10所示。

图9 Gfosub子程序建立模型仿真接触力Fig.9 The force simulated by the Gfosub model

图10 ADAMS自身接触建立模型仿真接触力Fig.10 The force simulated by the model of ADAMS’s contact

由于ADAMS的Parasoilds求解器不能处理复杂形状的接触问题,而自身默认的求解器又不适合处理微观接触,误差较大,其仿真的接触力扩大十几倍,造成其摩擦拖动力也相应扩大,无法准确反映轻载、高速下轴承的打滑率情况,而且接触力波动较大;而通过对ADAMS进行二次开发,利用Gfosub子程序建立的轴承动力学模型的仿真精度较高,接触力的最大值为141.59 N,利用相关经验公式进行预估[21],得到理论接触力为142.85 N,两者误差为0.82%,在可接受范围内,从而证明了子程序计算的优越性与可行性。

图11 保持架转速测量验证实验台Fig.11 A verification experiment platform of cage rotation speed

3 实验验证

3.1 转速测量实验台的搭建

为了验证所建立的轴承动力学模型的正确性,搭建了高速滚动轴承保持架转速测量实验台,如图11所示。整体结构由机械主轴-轴承系统和测试系统两部分构成,其中测试系统主要包括用于发射和接收红外信号的激光转速传感器、用于反射激光信号的反光条、用于采集脉冲信号的数据采集系统以及用于信号处理的计算机等。在测量之前,先将反光片贴在保持架的轴向端面上,再将激光速度传感器固定在磁铁支架,并与数据采集系统和计算机相连,调整传感器的位置,确保发出的光源能够直射到反光条上。在测量过程中,保持架每旋转一周,激光转速传感器即可接收一次反射的红外信号,进而产生一个电压脉冲,最终经过计算机进行相关数据处理便可得到保持架实际转速。

为了检验该装置测量转速的精度,先以已知转速的机械轴为测量对象预先做了一组实验。将机械轴由停机状态下匀加速到12 000 r/min,通过机械轴的实际转速和测量转速的对比,发现误差小于1%,说明该装置测量精度满足验证要求。

3.2 实验结果对比分析

基于上述二次开发方案,对7008C角接触球轴承动态特性开展仿真与实验对比验证,分析了轴向预紧力为50、100、150、200 N,内圈驱动转速在2 000、4 000、6 000 r/min工况下的保持架转速,结果如图12所示。

(a)转速2 000 r/min

(b)转速4 000 r/min

(c)转速6 000 r/min图12 保持架实际转速测量与仿真对比Fig.12 A comparison of actual measurement and simulation of cage speed

观察图12可以发现:本文建立的角接触球轴承动力学模型在2 000 r/min时,预紧力对保持架转速影响较弱,且保持架转速的仿真值与实验值较为接近,最大相差0.236 3 rad/s,此时相对误差为0.26%。在4 000 r/min下,预紧力对保持架转速的影响微小,随着预紧力的增加,保持架转速逐渐增加,而保持架转速的仿真值与实验值极为接近,此时最大相差0.227 0 rad/s,相对误差为0.125%。在6 000 r/min下,预紧力对保持架的转速影响较小,而保持架转速的仿真值与实验值最大相差0.630 3 rad/s,此时相对误差为0.23%。实验测量值与模型仿真值的相差不超过0.5%,且误差产生的原因可能是模型中些许参数存在相关误差,或保持架在实际测量过程中存在测量误差,但整体上证明建立的轴承模型运动分析的准确性。

4 结 论

针对ADAMS自身接触器计算精度不足的问题,本文基于ADAMS软件自身动力学分析优势,对其接触求解器进行了相关开发,建立了角接触球轴承动力学模型,对轴承之间的载荷分布与接触进行了精确建模,并利用实验验证了模型的仿真计算精度,具体结论如下:

(1)基于轴承全局坐标系,利用Gfosub子程序编写轴承零件之间的作用力子程序,根据求解方法的特点应用GSTIFF I3进行求解,实现了面向ADAMS软件二次开发的接触与运动精确建模与计算,经过与经验公式对比,接触力仿真误差不超过1%;

(2)通过搭建保持架转速测量实验台,对建立的动力学模型进行相关验证,发现在低速时,模型能够准确地仿真计算出保持架的转速,且仿真与实验测量误差不超过0.5%,证明了模型的可靠性;

(3)在高速时,建立的模型仿真得到的保持架转速与实验测得转速存在着一定的差别,还有待于进一步细化研究。