基于ABAQUS的关节轴承径向受载试验的热力耦合分析

霍亚军，尹忠慰，李虎林，曲艳峰

(上海交通大学 机械与动力工程学院，上海 200240)

关节轴承是一种球面滑动轴承，主要由外球面内圈和内球面外圈组成[1]。因为关节轴承的球形滑动接触面大，允许的倾斜角大，且多数关节轴承采用了表面磷化、镀铬、滑动面衬里及镶垫等特殊的工艺处理，所以具有结构简单、承受载荷大、抗冲击能力强、抗腐蚀、自调心、润滑好、耐磨损和长寿命等特点，广泛用于工程机械、工程结构、轻工机械、水利机械、军工机械等[1-4]。随着工业的飞速发展和新装备的研发，特别是在军事、航空航天和尖端科技领域，对关节轴承的性能提出了更高的要求。

目前，研究人员主要通过理论分析计算、试验和数值仿真的方法对关节轴承进行研究。文献[5]运用有限元法建立了关节轴承的仿真模型，并通过激光散射光弹法(SLP)测得环氧树脂模型的应力场，通过对比光弹性条纹与仿真结果，验证了有限元数值仿真模型的有效性，分析了配合间隙、轴承座是否为刚体等对试验和仿真结果的影响，并得出了经验公式假设的接触应力分布情况与实际不符的结论。文献[6]对自润滑关节轴承进行了三维顺序耦合热应力分析，并将仿真分析结果与相应的试验结果进行对比，验证了有限元分析方法的可行性，并分析了温度对接触应力的影响。文献[7]对大型推力关节轴承进行有限元仿真分析，研究了关节轴承在轴向载荷、纯径向载荷和轴、径向复合载荷3种情况下接触应力、球体内孔径向变形及轴向变形的变化规律，指出了内孔径的最大变形可作为确定轴承与轴配合间隙的理论依据。文献[8]利用ANSYS软件对关节轴承GEZ101ES进行结构优化，在装配尺寸不变的条件下，得出当球径尺寸减少1 mm时，最大压应力最小，并通过试验验证了优化后的关节轴承磨损寿命得到增长。文献[9]对自润滑向心关节轴承摩擦磨损寿命模型进行了分析研究，通过试验发现，关节轴承摩擦表面附近的温度随摆动频率的增大而不断升高，且3种不同运动形式对温度的影响也不同，其中复合摆动条件下的温度最大，旋转摆动次之，而倾斜摆动条件下温度值最低，并指出更全面分析出现这种差异的原因还需要从有限元热分析等方面进行考察。

关节轴承工作时内圈外球面与外圈内球面之间发生相对滑动，并通过该球面副传递载荷，所以内、外圈之间的接触应力对关节轴承的运行性能具有重要影响。此外，接触应力直接影响内外球面间的摩擦力，随着工作时间的增加，其对关节轴承的温升和接触面的磨损量(寿命)影响逐渐显现，并且二者的变化又会影响关节轴承工作面接触应力，所以关节轴承工作过程中，温度场和应力场是相互耦合的。但目前针对关节轴承的有限元分析主要是静力分析或者顺序耦合热力分析。针对上述问题，基于ABAQUS编写了用于热力耦合分析的用户子程序，对自润滑向心关节轴承进行完全耦合热力分析，对温度场和应力场同时求解。

1 理论基础与仿真分析过程

1.1 热力耦合分析方程[10-11]

对于连续介质，其能量守恒方程为

(1)

式中：Ω为体积；S为边界；ρ为介质密度；vi为速度场；U为单位质量介质内能；Q为单位体积热流；bi为单位体积力；Pi为单位面积上的边界力；H为边界热流密度。

根据动量守恒，建立介质的力平衡方程为

(2)

引入柯西应力分量σij(j=1,2,3,其值分别对应直角坐标系中的x,y,z方向)，压力用柯西应力表示为

Pi=niσij，

(3)

式中：ni为边界法线方向。

结合(1)～(3)式可得热力耦合的能量守恒方程为

(4)

根据虚功原理，可建立结构位移ui所满足的关系

(5)

式中：δui表示虚位移;xj表示方向轴。

单元位移矢量uE和节点位移矢量uN的关系为

uE(x,t)=N(x)uN(t)，

(6)

式中：N(x)为形函数矩阵；t为时间变量。上式对时间求导得形变率为

(7)

对于温度场，则有

TE(x,t)=B(x)TN(t)，

(8)

式中：TE(x,t)为单元温度矢量；B(x)为温度场的插值函数矩阵；TN(t)为节点温度矢量。

单元应变矢量为

ε(x,t)=LuE(x,t)=G(x)uN(t)，

(9)

式中：L为微分算子；G(x)为单元应变与节点位移之间关系的几何矩阵。

温度梯度矢量为

(10)

式中：A(x)为单元温度梯度矢量与节点温度矢量之间的关系矩阵。

结构瞬态温度场和热应力应变场分析的有限元方程为

(11)

KTTN(t)]=0，

(12)

式中：Ku，MT分别为结构刚度矩阵、热学刚度矩阵；F(t)为受力向量；Cu为热熔矩阵；KT为热传导矩阵；Mu为热力耦合矩阵；D为耗散向量；R为热载荷向量。合并(11)和(12)式可得

(13)

Z(t)=D+R+KTTN(t)。

对于接触摩擦生热问题，热载荷矢量R可表示为

R=κFfvr，

(14)

式中：κ为热功转换系数；Ff为摩擦力；vr为表面相对滑动速度。

1.2 ABAQUS完全耦合热力分析

热力耦合分析可分为顺序耦合热力分析和完全耦合热力分析。顺序耦合热力分析是首先进行传热分析，然后将得到的温度场作为边界条件，进行静力分析，得到应力应变场；完全耦合热力分析是考虑了温度场和应力应变场之间的耦合作用，对温度场和应力应变场同时进行求解。使用ABAQUS软件，对自润滑向心关节轴承进行完全耦合热力分析的主要过程如下：

(1)分析过程。通过前处理，得到有限元分析模型，包括几何模型的建立、材料参数的定义、装配并划分网格、分析步的设定、接触关系的设定、边界载荷的施加；然后建立分析作业，提交后进行相关计算，最后通过后处理模块(Visualiztion)查看分析结果。自润滑向心关节轴承的几何结构与有限元模型如图1所示[12]。

图1 关节轴承的几何结构与有限元模型

在建立有限元模型的过程中，忽略对计算结果影响不大的轴承实体模型的细微结构，如倒角、圆角等。对于衬垫层，将有限元模型的外圈分割成2部分，分别设置衬垫和外圈的材料属性[6]。内、外圈材料为铝合金，衬垫基体材料为PTFE纤维与聚酰胺纤维编织物[13]，衬垫看作是宏观各向同性材料[6]，材料的力学性能见表1[14-16]。通常材料的热学性能参数(热膨胀系数α、热传导系数λ、比热C)会随着温度的变化而变化，但在分析过程中进行简化处理，将其设定为常数。下标1表示沿衬垫厚度方向，2，3表示在衬垫平面内。

表1 关节轴承的材料参数

采用扫掠网格划分技术，选择线性减缩积分单元C3D8RT，建立Coupled Temp-Displacement分析步，在Response选项中选择Steady-State计算自润滑关节轴承的稳态温度场分布[17]。接触关系属性中需要设定接触面的热传导系数和摩擦生成的热量在主从接触面的分配情况[18]。此外，在接触关系定义中，将表面散热系数设定为关节轴承内、外圈材料与关节轴承夹具材料热传导系数的均值[19]。接触关系属性参数同样也会随着温度的变化而变化，文中将其设为恒定值。分别将内圈内圆柱面和外圈外圆柱面设定为刚性体，通过参考点控制其自由度，进而施加力和位移边界条件。关节轴承受到的径向载荷为510 kN，关节轴承的摆动角度为±25°、摆动频率为每分钟10次。通过编写用户子程序实现完全耦合热力分析，子程序中采用(14)式计算摩擦能耗，关节轴承的运动参数可直接代入该式进行计算，故在施加位移边界条件时可固定约束内圈绕轴线摆动的自由度。

(2)用户子程序。ABAQUS提供了大量的用户子程序(User Subroutines)作为二次开发的平台，用户可根据自己的需要定义符合特定问题的模型[20]。通常，在ABAQUS中使用ABAQUS/Standard求解器进行热力耦合分析时，如果选择稳态分析(Steady-State)，求解器会自动忽略所有节点的温度自由度，但可以通过编写用户子程序FRIC来计算每一增量步的摩擦能耗SPD，从而计算产生的热量。开发用于完全耦合热力分析的子程序有2种方法：(1)通过编写子程序FRIC实现某类型摩擦定律，如库仑摩擦定律，此时可以定义切向应力(摩擦应力TAU)来计算摩擦能耗SPD，即摩擦模型法；(2)根据(14)式计算摩擦能耗SPD，即公式法。采用公式法时，可以通过子程序中传递的节点接触正应力CPRESS和每个增量步时间内该节点滑动的位移来计算摩擦能耗SPD，热功转换系数κ取0.85。文中采用公式法进行分析，用户子程序FRIC的详细介绍可参考文献[21]。

2 结果分析与讨论

自润滑向心关节轴承在不同初始温度条件下的稳态温度场如图2所示。图2a和图2c分别为环境温度18 ℃和7 ℃时轴承整体稳态温度场，其温度最大值分别为107.7 ℃和96.73 ℃，最大温度均处于内圈与衬垫的接触面上;图2b和图2d分别为环境温度18 ℃和7 ℃时外圈稳态温度场，温度最大值分别为69.99 ℃和58.99 ℃，最大温度均处于外圈与衬垫的接触面上。从图2还可以看出，轴承温度分布范围较大，接触中心处温度较高，向两侧逐渐降低，这主要受关节轴承球面副接触应力分布情况及仿真分析中衬垫层与外圈接触关系设定的影响。

图2 关节轴承稳态温度场

通常试验使用接触式温度传感器测量内圈或外圈端面附近温度[22-23]，通过把相应节点温度值与试验测量的温度值进行对比，可以说明仿真计算结果的有效性。在与仿真分析相对应的试验中，用铂电阻通过XMZ数字显示仪表测量轴承外圈的端面温度，精度误差不超过±1%，可测量轴承端面温度范围为0～300 ℃[22]。与试验测量值[12]和顺序耦合热力分析仿真值[6]的对比情况见表2。由表可知，完全耦合热力分析得到的温度值更接近试验测量值，相对误差为3.3%，小于顺序耦合热力分析的相对误差(8.4%)。所以，基于ABAQUS通过二次开发热力耦合分析子程序，可以有效地对自润滑向心关节轴承进行完全耦合热力分析。仿真值小于试验测量值的原因可能是：对衬垫的热学性能参数进行了简化；材料参数和某些边界条件设定为常值，忽略了温度对其的影响；没有考虑衬垫层磨损的影响。

表2 分析结果对比

环境温度为18 ℃时，自润滑向心关节轴承内、外圈和衬垫的Mises应力如图3所示。内、外圈的最大Mises应力值分别为356.2和287 MPa，分别位于内、外径表面上；衬垫的最大Mises应力为135.7 MPa，位于衬垫与内圈接触面上。

图3 关节轴承热应力分布图

环境温度为18 ℃时，关节轴承沿径向力方向(x轴)的位移分布如图4所示。由图4a和图4b可知，轴承沿x轴方向的最大位移为-0.033 53 mm，位于内圈接触面上。由图4c可知，在不考虑热影响的情况下时，内圈接触面上的位移为-0.038 21 mm，绝对值大于0.033 53 mm，说明轴承温升引起的热膨胀抵消了部分弹性变形。

图4 关节轴承沿径向力方向的位移分布图

环境温度为18 ℃时，关节轴承接触应力分布如图5所示。其中图5a为温度达到稳态后接触应力分布情况，最大接触应力为214.6 MPa；为了分析温度因素对接触应力的影响，建立了如图5b所示的内圈外球面上的节点路径Path-1，该节点路径沿接触的圆周方向分布，包角为180°。

图5 关节轴承接触应力分布图

沿路径Path-1完全耦合热力分析和顺序耦合热力分析的接触应力变化情况如图6所示。由图可知，同一位置处完全耦合热力分析的接触应力大于顺序耦合热力分析的接触应力，随着沿z向位置增大，接触应力均呈先增大后减小的趋势。在远离轴承接触中心两侧(z=0和z=126 mm处)，接触应力增加最多，增加了64.66 MPa，增幅为180%；随着位置向接触中心(z=63 mm处)移动，接触应力增加值逐渐减小，最大接触应力增加了39.20 MPa，增幅为22%。由此可知，温度对关节轴承的位移场与接触应力影响较大。

图6 接触应力沿路径Path-1的分布情况(环境温度18 ℃)

3 结论

(1)通过完全耦合热力分析得到了关节轴承的稳态温度场，其轴承外圈的最大温度值较顺序耦合热力分析更接近于试验测量值。

(2)与顺序耦合热力分析相比，温度变化引起的热膨胀抵消了部分弹性变形量，使得内圈沿径向力方向的最大位移有所减小，而接触应力和最大接触应力均有所增加。