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

(上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院，上海 200240)

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

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

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

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

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

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

(1)

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

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

(2)

引入柯西應(yīng)力分量σij(j=1,2,3,其值分別對(duì)應(yīng)直角坐標(biāo)系中的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é)點(diǎn)位移矢量uN的關(guān)系為

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

(6)

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

(7)

對(duì)于溫度場(chǎng)，則有

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

(8)

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

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

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

(9)

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

溫度梯度矢量為

(10)

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

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

(11)

KTTN(t)]=0，

(12)

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

(13)

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

對(duì)于接觸摩擦生熱問(wèn)題，熱載荷矢量R可表示為

R=κFfvr，

(14)

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

1.2 ABAQUS完全耦合熱力分析

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

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

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

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

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

采用掃掠網(wǎng)格劃分技術(shù)，選擇線性減縮積分單元C3D8RT，建立Coupled Temp-Displacement分析步，在Response選項(xiàng)中選擇Steady-State計(jì)算自潤(rùn)滑關(guān)節(jié)軸承的穩(wěn)態(tài)溫度場(chǎng)分布[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ù)同樣也會(huì)隨著溫度的變化而變化，文中將其設(shè)為恒定值。分別將內(nèi)圈內(nèi)圓柱面和外圈外圓柱面設(shè)定為剛性體，通過(guò)參考點(diǎn)控制其自由度，進(jìn)而施加力和位移邊界條件。關(guān)節(jié)軸承受到的徑向載荷為510 kN，關(guān)節(jié)軸承的擺動(dòng)角度為±25°、擺動(dòng)頻率為每分鐘10次。通過(guò)編寫用戶子程序?qū)崿F(xiàn)完全耦合熱力分析，子程序中采用(14)式計(jì)算摩擦能耗，關(guān)節(jié)軸承的運(yùn)動(dòng)參數(shù)可直接代入該式進(jìn)行計(jì)算，故在施加位移邊界條件時(shí)可固定約束內(nèi)圈繞軸線擺動(dòng)的自由度。

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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