霍亞軍，尹忠慰，李虎林，曲艷峰

(上海交通大學 機械與動力工程學院，上海 200240)

關(guān)節(jié)軸承是一種球面滑動軸承，主要由外球面內(nèi)圈和內(nèi)球面外圈組成[1]。因為關(guān)節(jié)軸承的球形滑動接觸面大，允許的傾斜角大，且多數(shù)關(guān)節(jié)軸承采用了表面磷化、鍍鉻、滑動面襯里及鑲墊等特殊的工藝處理，所以具有結(jié)構(gòu)簡單、承受載荷大、抗沖擊能力強、抗腐蝕、自調(diào)心、潤滑好、耐磨損和長壽命等特點，廣泛用于工程機械、工程結(jié)構(gòu)、輕工機械、水利機械、軍工機械等[1-4]。隨著工業(yè)的飛速發(fā)展和新裝備的研發(fā)，特別是在軍事、航空航天和尖端科技領(lǐng)域，對關(guān)節(jié)軸承的性能提出了更高的要求。

目前，研究人員主要通過理論分析計算、試驗和數(shù)值仿真的方法對關(guān)節(jié)軸承進行研究。文獻[5]運用有限元法建立了關(guān)節(jié)軸承的仿真模型，并通過激光散射光彈法(SLP)測得環(huán)氧樹脂模型的應(yīng)力場，通過對比光彈性條紋與仿真結(jié)果，驗證了有限元數(shù)值仿真模型的有效性，分析了配合間隙、軸承座是否為剛體等對試驗和仿真結(jié)果的影響，并得出了經(jīng)驗公式假設(shè)的接觸應(yīng)力分布情況與實際不符的結(jié)論。文獻[6]對自潤滑關(guān)節(jié)軸承進行了三維順序耦合熱應(yīng)力分析，并將仿真分析結(jié)果與相應(yīng)的試驗結(jié)果進行對比，驗證了有限元分析方法的可行性，并分析了溫度對接觸應(yīng)力的影響。文獻[7]對大型推力關(guān)節(jié)軸承進行有限元仿真分析，研究了關(guān)節(jié)軸承在軸向載荷、純徑向載荷和軸、徑向復(fù)合載荷3種情況下接觸應(yīng)力、球體內(nèi)孔徑向變形及軸向變形的變化規(guī)律，指出了內(nèi)孔徑的最大變形可作為確定軸承與軸配合間隙的理論依據(jù)。文獻[8]利用ANSYS軟件對關(guān)節(jié)軸承GEZ101ES進行結(jié)構(gòu)優(yōu)化，在裝配尺寸不變的條件下，得出當球徑尺寸減少1 mm時，最大壓應(yīng)力最小，并通過試驗驗證了優(yōu)化后的關(guān)節(jié)軸承磨損壽命得到增長。文獻[9]對自潤滑向心關(guān)節(jié)軸承摩擦磨損壽命模型進行了分析研究，通過試驗發(fā)現(xiàn)，關(guān)節(jié)軸承摩擦表面附近的溫度隨擺動頻率的增大而不斷升高，且3種不同運動形式對溫度的影響也不同，其中復(fù)合擺動條件下的溫度最大，旋轉(zhuǎn)擺動次之，而傾斜擺動條件下溫度值最低，并指出更全面分析出現(xiàn)這種差異的原因還需要從有限元熱分析等方面進行考察。

關(guān)節(jié)軸承工作時內(nèi)圈外球面與外圈內(nèi)球面之間發(fā)生相對滑動，并通過該球面副傳遞載荷，所以內(nèi)、外圈之間的接觸應(yīng)力對關(guān)節(jié)軸承的運行性能具有重要影響。此外，接觸應(yīng)力直接影響內(nèi)外球面間的摩擦力，隨著工作時間的增加，其對關(guān)節(jié)軸承的溫升和接觸面的磨損量(壽命)影響逐漸顯現(xiàn)，并且二者的變化又會影響關(guān)節(jié)軸承工作面接觸應(yīng)力，所以關(guān)節(jié)軸承工作過程中，溫度場和應(yīng)力場是相互耦合的。但目前針對關(guān)節(jié)軸承的有限元分析主要是靜力分析或者順序耦合熱力分析。針對上述問題，基于ABAQUS編寫了用于熱力耦合分析的用戶子程序，對自潤滑向心關(guān)節(jié)軸承進行完全耦合熱力分析，對溫度場和應(yīng)力場同時求解。

1 理論基礎(chǔ)與仿真分析過程

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

對于連續(xù)介質(zhì)，其能量守恒方程為

(1)

式中：Ω為體積；S為邊界；ρ為介質(zhì)密度；vi為速度場；U為單位質(zhì)量介質(zhì)內(nèi)能；Q為單位體積熱流；bi為單位體積力；Pi為單位面積上的邊界力；H為邊界熱流密度。

根據(jù)動量守恒，建立介質(zhì)的力平衡方程為

(2)

引入柯西應(yīng)力分量σij(j=1,2,3,其值分別對應(yīng)直角坐標系中的x,y,z方向)，壓力用柯西應(yīng)力表示為

Pi=niσij，

(3)

式中：ni為邊界法線方向。

結(jié)合(1)～(3)式可得熱力耦合的能量守恒方程為

(4)

根據(jù)虛功原理，可建立結(jié)構(gòu)位移ui所滿足的關(guān)系

(5)

式中：δui表示虛位移;xj表示方向軸。

單元位移矢量uE和節(jié)點位移矢量uN的關(guān)系為

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

(6)

式中：N(x)為形函數(shù)矩陣；t為時間變量。上式對時間求導(dǎo)得形變率為

(7)

對于溫度場，則有

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

(8)

式中：TE(x,t)為單元溫度矢量；B(x)為溫度場的插值函數(shù)矩陣；TN(t)為節(jié)點溫度矢量。

單元應(yīng)變矢量為

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

(9)

式中：L為微分算子；G(x)為單元應(yīng)變與節(jié)點位移之間關(guān)系的幾何矩陣。

溫度梯度矢量為

(10)

式中：A(x)為單元溫度梯度矢量與節(jié)點溫度矢量之間的關(guān)系矩陣。

結(jié)構(gòu)瞬態(tài)溫度場和熱應(yīng)力應(yīng)變場分析的有限元方程為

(11)

KTTN(t)]=0，

(12)

式中：Ku，MT分別為結(jié)構(gòu)剛度矩陣、熱學剛度矩陣；F(t)為受力向量；Cu為熱熔矩陣；KT為熱傳導(dǎo)矩陣；Mu為熱力耦合矩陣；D為耗散向量；R為熱載荷向量。合并(11)和(12)式可得

(13)

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

對于接觸摩擦生熱問題，熱載荷矢量R可表示為

R=κFfvr，

(14)

式中：κ為熱功轉(zhuǎn)換系數(shù)；Ff為摩擦力；vr為表面相對滑動速度。

1.2 ABAQUS完全耦合熱力分析

熱力耦合分析可分為順序耦合熱力分析和完全耦合熱力分析。順序耦合熱力分析是首先進行傳熱分析，然后將得到的溫度場作為邊界條件，進行靜力分析，得到應(yīng)力應(yīng)變場；完全耦合熱力分析是考慮了溫度場和應(yīng)力應(yīng)變場之間的耦合作用，對溫度場和應(yīng)力應(yīng)變場同時進行求解。使用ABAQUS軟件，對自潤滑向心關(guān)節(jié)軸承進行完全耦合熱力分析的主要過程如下：

(1)分析過程。通過前處理，得到有限元分析模型，包括幾何模型的建立、材料參數(shù)的定義、裝配并劃分網(wǎng)格、分析步的設(shè)定、接觸關(guān)系的設(shè)定、邊界載荷的施加；然后建立分析作業(yè)，提交后進行相關(guān)計算，最后通過后處理模塊(Visualiztion)查看分析結(jié)果。自潤滑向心關(guān)節(jié)軸承的幾何結(jié)構(gòu)與有限元模型如圖1所示[12]。

圖1 關(guān)節(jié)軸承的幾何結(jié)構(gòu)與有限元模型

在建立有限元模型的過程中，忽略對計算結(jié)果影響不大的軸承實體模型的細微結(jié)構(gòu)，如倒角、圓角等。對于襯墊層，將有限元模型的外圈分割成2部分，分別設(shè)置襯墊和外圈的材料屬性[6]。內(nèi)、外圈材料為鋁合金，襯墊基體材料為PTFE纖維與聚酰胺纖維編織物[13]，襯墊看作是宏觀各向同性材料[6]，材料的力學性能見表1[14-16]。通常材料的熱學性能參數(shù)(熱膨脹系數(shù)α、熱傳導(dǎo)系數(shù)λ、比熱C)會隨著溫度的變化而變化，但在分析過程中進行簡化處理，將其設(shè)定為常數(shù)。下標1表示沿襯墊厚度方向，2，3表示在襯墊平面內(nèi)。

表1 關(guān)節(jié)軸承的材料參數(shù)

采用掃掠網(wǎng)格劃分技術(shù)，選擇線性減縮積分單元C3D8RT，建立Coupled Temp-Displacement分析步，在Response選項中選擇Steady-State計算自潤滑關(guān)節(jié)軸承的穩(wěn)態(tài)溫度場分布[17]。接觸關(guān)系屬性中需要設(shè)定接觸面的熱傳導(dǎo)系數(shù)和摩擦生成的熱量在主從接觸面的分配情況[18]。此外，在接觸關(guān)系定義中，將表面散熱系數(shù)設(shè)定為關(guān)節(jié)軸承內(nèi)、外圈材料與關(guān)節(jié)軸承夾具材料熱傳導(dǎo)系數(shù)的均值[19]。接觸關(guān)系屬性參數(shù)同樣也會隨著溫度的變化而變化，文中將其設(shè)為恒定值。分別將內(nèi)圈內(nèi)圓柱面和外圈外圓柱面設(shè)定為剛性體，通過參考點控制其自由度，進而施加力和位移邊界條件。關(guān)節(jié)軸承受到的徑向載荷為510 kN，關(guān)節(jié)軸承的擺動角度為±25°、擺動頻率為每分鐘10次。通過編寫用戶子程序?qū)崿F(xiàn)完全耦合熱力分析，子程序中采用(14)式計算摩擦能耗，關(guān)節(jié)軸承的運動參數(shù)可直接代入該式進行計算，故在施加位移邊界條件時可固定約束內(nèi)圈繞軸線擺動的自由度。

(2)用戶子程序。ABAQUS提供了大量的用戶子程序(User Subroutines)作為二次開發(fā)的平臺，用戶可根據(jù)自己的需要定義符合特定問題的模型[20]。通常，在ABAQUS中使用ABAQUS/Standard求解器進行熱力耦合分析時，如果選擇穩(wěn)態(tài)分析(Steady-State)，求解器會自動忽略所有節(jié)點的溫度自由度，但可以通過編寫用戶子程序FRIC來計算每一增量步的摩擦能耗SPD，從而計算產(chǎn)生的熱量。開發(fā)用于完全耦合熱力分析的子程序有2種方法：(1)通過編寫子程序FRIC實現(xiàn)某類型摩擦定律，如庫侖摩擦定律，此時可以定義切向應(yīng)力(摩擦應(yīng)力TAU)來計算摩擦能耗SPD，即摩擦模型法；(2)根據(jù)(14)式計算摩擦能耗SPD，即公式法。采用公式法時，可以通過子程序中傳遞的節(jié)點接觸正應(yīng)力CPRESS和每個增量步時間內(nèi)該節(jié)點滑動的位移來計算摩擦能耗SPD，熱功轉(zhuǎn)換系數(shù)κ取0.85。文中采用公式法進行分析，用戶子程序FRIC的詳細介紹可參考文獻[21]。

2 結(jié)果分析與討論

自潤滑向心關(guān)節(jié)軸承在不同初始溫度條件下的穩(wěn)態(tài)溫度場如圖2所示。圖2a和圖2c分別為環(huán)境溫度18 ℃和7 ℃時軸承整體穩(wěn)態(tài)溫度場，其溫度最大值分別為107.7 ℃和96.73 ℃，最大溫度均處于內(nèi)圈與襯墊的接觸面上;圖2b和圖2d分別為環(huán)境溫度18 ℃和7 ℃時外圈穩(wěn)態(tài)溫度場，溫度最大值分別為69.99 ℃和58.99 ℃，最大溫度均處于外圈與襯墊的接觸面上。從圖2還可以看出，軸承溫度分布范圍較大，接觸中心處溫度較高，向兩側(cè)逐漸降低，這主要受關(guān)節(jié)軸承球面副接觸應(yīng)力分布情況及仿真分析中襯墊層與外圈接觸關(guān)系設(shè)定的影響。

圖2 關(guān)節(jié)軸承穩(wěn)態(tài)溫度場

通常試驗使用接觸式溫度傳感器測量內(nèi)圈或外圈端面附近溫度[22-23]，通過把相應(yīng)節(jié)點溫度值與試驗測量的溫度值進行對比，可以說明仿真計算結(jié)果的有效性。在與仿真分析相對應(yīng)的試驗中，用鉑電阻通過XMZ數(shù)字顯示儀表測量軸承外圈的端面溫度，精度誤差不超過±1%，可測量軸承端面溫度范圍為0～300 ℃[22]。與試驗測量值[12]和順序耦合熱力分析仿真值[6]的對比情況見表2。由表可知，完全耦合熱力分析得到的溫度值更接近試驗測量值，相對誤差為3.3%，小于順序耦合熱力分析的相對誤差(8.4%)。所以，基于ABAQUS通過二次開發(fā)熱力耦合分析子程序，可以有效地對自潤滑向心關(guān)節(jié)軸承進行完全耦合熱力分析。仿真值小于試驗測量值的原因可能是：對襯墊的熱學性能參數(shù)進行了簡化；材料參數(shù)和某些邊界條件設(shè)定為常值，忽略了溫度對其的影響；沒有考慮襯墊層磨損的影響。

表2 分析結(jié)果對比

環(huán)境溫度為18 ℃時，自潤滑向心關(guān)節(jié)軸承內(nèi)、外圈和襯墊的Mises應(yīng)力如圖3所示。內(nèi)、外圈的最大Mises應(yīng)力值分別為356.2和287 MPa，分別位于內(nèi)、外徑表面上；襯墊的最大Mises應(yīng)力為135.7 MPa，位于襯墊與內(nèi)圈接觸面上。

圖3 關(guān)節(jié)軸承熱應(yīng)力分布圖

環(huán)境溫度為18 ℃時，關(guān)節(jié)軸承沿徑向力方向(x軸)的位移分布如圖4所示。由圖4a和圖4b可知，軸承沿x軸方向的最大位移為-0.033 53 mm，位于內(nèi)圈接觸面上。由圖4c可知，在不考慮熱影響的情況下時，內(nèi)圈接觸面上的位移為-0.038 21 mm，絕對值大于0.033 53 mm，說明軸承溫升引起的熱膨脹抵消了部分彈性變形。

圖4 關(guān)節(jié)軸承沿徑向力方向的位移分布圖

環(huán)境溫度為18 ℃時，關(guān)節(jié)軸承接觸應(yīng)力分布如圖5所示。其中圖5a為溫度達到穩(wěn)態(tài)后接觸應(yīng)力分布情況，最大接觸應(yīng)力為214.6 MPa；為了分析溫度因素對接觸應(yīng)力的影響，建立了如圖5b所示的內(nèi)圈外球面上的節(jié)點路徑Path-1，該節(jié)點路徑沿接觸的圓周方向分布，包角為180°。

圖5 關(guān)節(jié)軸承接觸應(yīng)力分布圖

沿路徑Path-1完全耦合熱力分析和順序耦合熱力分析的接觸應(yīng)力變化情況如圖6所示。由圖可知，同一位置處完全耦合熱力分析的接觸應(yīng)力大于順序耦合熱力分析的接觸應(yīng)力，隨著沿z向位置增大，接觸應(yīng)力均呈先增大后減小的趨勢。在遠離軸承接觸中心兩側(cè)(z=0和z=126 mm處)，接觸應(yīng)力增加最多，增加了64.66 MPa，增幅為180%；隨著位置向接觸中心(z=63 mm處)移動，接觸應(yīng)力增加值逐漸減小，最大接觸應(yīng)力增加了39.20 MPa，增幅為22%。由此可知，溫度對關(guān)節(jié)軸承的位移場與接觸應(yīng)力影響較大。

圖6 接觸應(yīng)力沿路徑Path-1的分布情況(環(huán)境溫度18 ℃)

3 結(jié)論

(1)通過完全耦合熱力分析得到了關(guān)節(jié)軸承的穩(wěn)態(tài)溫度場，其軸承外圈的最大溫度值較順序耦合熱力分析更接近于試驗測量值。

(2)與順序耦合熱力分析相比，溫度變化引起的熱膨脹抵消了部分彈性變形量，使得內(nèi)圈沿徑向力方向的最大位移有所減小，而接觸應(yīng)力和最大接觸應(yīng)力均有所增加。