深沟球轴承运转过程动态特性有限元分析_林腾蛟

振　动　与　冲　击第28卷第1期J O U R N A LO FV I B R A T I O NA N DS H O C KV o l .28N o .12009　深沟球轴承运转过程动态特性有限元分析基金项目:国家自然科学基金资助项目(编号:50675232,50735008);长江学者和创新团队发展计划资助项目(编号:I R T 0763)收稿日期:2007-11-12　修改稿收到日期:2008-03-20第一作者林腾蛟男,博士后,教授,1968年5月生林腾蛟,荣　崎,李润方,邵毅敏(重庆大学机械传动国家重点实验室,重庆　400044) 摘　要:综合考虑轴承径向载荷及转速的影响,应用A N S Y S /L S -D Y N A 软件建立了深沟球轴承多体动力接触有限元模型;以显式动力学有限元法为基础,采用全积分单元算法控制沙漏,设置质量缩放系数缩减计算时间,对内圈施加不同转速时的深沟球轴承进行动力接触分析,得出了深沟球轴承运转过程的动态响应及滚动体的应力分布。

滚动体的最大和最小线速度分别出现在与内、外圈接触点上,各转速下滚动体与内圈接触的应力基本相同,滚动体与外圈接触的应力随转速增高而相应增大,滚动体与外圈间接触力的波动大于内圈,而滚动体与保持架间的作用力较小。

研究表明,A N -S Y S /L S -D Y N A 是分析轴承运转过程动力接触问题十分有效的工具。

关键词:球轴承;动力接触;显式算法;有限元分析中图分类号:T H 113.1 文献标识码:A 滚动轴承具有摩擦力矩低、功率损耗小、维护方便、寿命长等特点,是应用极为广泛的重要机械基础件。

随着传动系统向高可靠和高速静音方向发展,对轴承的动力接触特性及振动噪声提出了愈来愈高的要求。

对滚动轴承的动力学分析始于19世纪70年代初。

W a l t e r 首先提出动力学设计分析方法,考虑轴承整个运动过程,建立动力学方程组,利用积分公式求得任一瞬间滚动体和保持架的位置、转速及轴承的内部滑动等。

而后G u p t a 进一步作了系统研究,但他们提出的动力学设计方法难以取得准确的分析解,必须应用计算机进行数值计算[1]。

随着计算机技术的快速发展,有限元仿真分析已广泛应用于解决工程问题。

目前国内外学者对滚动轴承的静力有限元分析较多,Y u a n K a n g 等借助接触有限元分析,修正了J o n e s 和H a r r i s 提出的深沟球轴承载荷与变形的关系,并通过了试验验证[2]。

唐云冰等建立了滚动轴承载荷分布的有限元分析模型,分析了载荷参数对轴承接触应力、接触角和变形的影响规律[3]。

刘宁等应用A N S Y S 对轧机滚动轴承进行接触应力分析,得出内圈与滚子之间接触应力的分布情况[4]。

而在滚动轴承动力学仿真方面,研究成果相对较少。

李松生等利用数值方法对超高速电主轴轴承的内部动力学状态进行计算机模拟仿真,分析了各种参数对轴承支承刚度的影响[5]。

樊莉等在A N S Y S /L S -D Y N A 中对滚动轴承的运动过程进行了数值模拟,讨论了滚动体接触应力的变化规律[6]。

笔者以深沟球轴承6206为例,综合考虑轴承径向载荷及转速的影响,将滚动体定义为线弹性材料模型,内、外圈及保持架定义为刚性材料模型,在A N S Y S 中建立深沟球轴承动力接触有限元模型,并采用显式动力学有限元法对深沟球轴承进行动力接触分析,求解轴承的位移、速度、加速度及滚动体单元应力。

1　深沟球轴承动力接触有限元模型1.1　轴承几何尺寸 表1给出了计算模型采用的轴承尺寸参数[7]。

表1　轴承参数符号描述值B 轴承宽度16m m D 球直径9.525m m D i 轴承内径30m m D o 轴承外径62m m d i 内滚道直径36.48m m d m 轴承平均直径46m m d o 外滚道直径55.53m m r i 内滚道曲率半径4.905m m r e 外滚道曲率半径4.953m m α接触角0°Z球数目91.2　单元类型选择及网格划分 建立深沟球轴承动力分析模型时,选用S O L I D 164单元,采用自由网格和映射网格相结合的方法进行网格划分,滚动体和内外圈采用六面体单元,保持架采用四面体单元,并将单元设置为全积分单元算法,进行沙漏控制。

由于A N S Y S /L S -D Y N A 中,S O L I D 164单元不具有旋转自由度,将轴承内圈表面定义为S H E L L 163单元,以便施加转速和载荷进行动力分析。

轴承有限元图1　轴承有限元网格网格如图1所示,共计有264935个节点,283035个单元。

1.3　材料模型选择A N S Y S 中定义刚体材料模型,可以大大减少计算时间,这是因为刚体内的所有节点的自由度都耦合到刚体的质心上,不论刚体定义了多少节点,刚体仅有6个自由度。

所有刚体的质量、质心和惯性都由刚性体体积和材料的密度自动计算得到,作用在刚性体上的力和力矩则由每个时间步的节点力和力矩合成,然后由此来计算刚性体的运动,节点位移由刚体质心的位移和旋转角度转换得到。

由于本文重点研究滚动体的动力学特性,在选择材料模型时,定义滚动体为线弹性材料模型,内、外圈及保持架为刚性材料模型。

内、外圈及滚动体的材料均用G C r 15钢,保持架为冷轧钢板,ρ=7830k g /m 3,E=2.07×1011P a ,ν=0.3。

1.4　边界条件设置轴承外圈固定,各刚性材料元件的约束形式如表2所示。

滚动体与内、外圈及保持架建立三组接触对,接触类型为自动面面接触,取静摩擦系数分别为0.1、0.1和0.1,动摩擦系数分别为0.002、0.002和0.002。

内圈施加径向载荷和转速,其中径向载荷F r 为5000N ,转速n i 分别为1000r /m i n 、2000r /m i n 、3000r /m i n 。

对滚动轴承同时施加径向载荷以及转速时,计算过程初始阶段通常不稳定。

为模拟轴承正常工作条件下的动态特性,先施加径向载荷,0.003s 时载荷引起的动态响应基本稳定后,再施加转速。

表2　刚体元件的约束情况平动约束转动约束内圈Z 方向X 、Y 方向外圈X 、Y 、Z 方向X 、Y 、Z 方向保持架X 、Y 、Z 方向X 、Y 方向2　L S -D Y N A 三维显式动力学算法 对于滚动轴承静态问题分析往往采用隐式积分算法,其优点是具有无条件稳定性,可取较大的时间步长。

然而,轴承运转过程仿真属于非线性多体动力接触问题,选择显式积分算法更适合于求解该类问题,主要体现在:1)滚动轴承接触物体多,计算模型复杂,动力仿真计算量大,隐式算法耗费机时多,对硬件要求也较高。

采用显式积分算法具有较高的计算效率,特别是可以通过合理设置质量缩放系数以减少计算时间。

2)L S -D Y N A 显式算法以中心差分法为主,不形成总体刚度矩阵,弹性项放在内力中,当质量矩阵和阻尼矩阵为对角阵时,不需要矩阵求逆,这对非线性分析非常有意义。

3)显式中心差分算法采用变步长积分法,每一时刻的积分步长由当前构形的稳定性条件控制,计算过程较为稳定,动力接触瞬时响应不连续性对计算结果的影响较小。

2.1　动力有限元基本方程 根据最小势能原理,运动微分方程的弱形式为[8]　∫V (ρu ··i -f i )δu i d V+∫V σi j δεi j d V -∫S σT i δu id S=0(1)式中:ρ为质量密度;u i 为节点位移张量;f i 为体力张量;δu i 为虚位移张量;σi j 为节点应力张量;δεi j 为节点虚应变应变张量;S σ为外力边界;T i 为面力张量。

将整个结构划分为一系列离散单元,用各单元的势能变分之和来近似的总势能,则∑eδu e T[∫VeρN TN d V u ··e+∫V eB Tσd V -∫V eN Tf d V-∫S σN TT d S ]=0(2)式中:N 为形函数;B 为应变矩阵;σ为应力向量;f 为体力向量;T 为面力向量;V e 为单元体积。

改写上式可得离散化运动方程为M u ··=P-F(3)其中:M 为总体质量矩阵M =∑e ∫Ve ρN TN d V (4)F 为应力散度向量,由单元应力场的等效节点力矢量组集而成F =∑e∫VeB Tσd V (5)P 为总体节点外载荷向量P =∑e(∫VeN T f d V+∫S σN TT d S )(6)L S -D Y N A 采用沙漏粘性阻力的办法来解决沙漏问题。

即在体单元各节点k 处,沿x i轴方向引入一个与沙漏模态变形方向相反的沙漏粘性阻尼力f i k =-αk ∑4j =1h i j Γj k ,i =1,2,3,k=1,2,…,8(7)式中:负号表示阻尼力与沙漏模态的变形方向相反;h i j =∑8k =1x ki Γj k 为沙漏模态的模;Γj k为第j 个沙漏基矢量的第k 个分量;系数αk =Q h g ρV 2/3e c /4,其中Q h g为用户定义常数,c 为材料声速。

将各节点的沙漏粘性阻尼力集成可得结构沙漏粘性阻尼力向量H ,同时考虑阻尼影响(阻尼矩阵C ),并119第1期 林腾蛟等:深沟球轴承运转过程动态特性有限元分析引入接触力R ,则动力接触问题离散化运动方程为M u ··=P-F+R+H -C u ·(8)2.2　显式算法与时步控制 时间积分采用显式中心差分方法,基本格式为u ··nt =M -1[P nt -F n t +R n t +H n t-C u ·n -12t ](9)u ·n +12t =u ·n -12t+u ··nt Δt n(10)u n +1t=u nt +u ·n +12t Δt n +12(11)式中:tn -12=(t n +t n -1)/2,tn +12=(t n +t n +1)/2,Δt n=(t n +1-t n -1)/2,Δt n +12=t n +1-t n;u ·n t ,u ··n +12t 和u n +1t 依次为t n时刻节点加速度向量,t n +12时刻的节点速度向量,t n +1时刻的节点位置坐标向量。

L S -D Y N A 采用变步长积分法,每一时刻的积分步长必须小于某个临界值,该值由网格中的最小单元决定,即Δt =m i n {Δt e 1,…,Δt e N }(12)式中:Δt e i为第i 个单元极限时步长;N 为单元总数。