基于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)與順序耦合熱力分析相比，溫度變化引起的熱膨脹抵消了部分彈性變形量，使得內圈沿徑向力方向的最大位移有所減小，而接觸應力和最大接觸應力均有所增加。