方群，李錦棒,，劉凱，王國(guó)榮，崔玉國(guó)

(1.寧波大學(xué) 機(jī)械工程與力學(xué)學(xué)院，浙江 寧波 315211；2.寧波達(dá)爾機(jī)械科技有限公司，浙江 寧波 315202；

3.西南石油大學(xué) 機(jī)電工程學(xué)院，成都 610500)

隨著機(jī)械設(shè)備向高速重載方向發(fā)展，服役軸承面臨著重大挑戰(zhàn)[1-4]。軸承在高速工況下的實(shí)際壽命往往遠(yuǎn)小于設(shè)計(jì)壽命，其接觸表面通常會(huì)出現(xiàn)早期的剝落失效[5-8]。研究發(fā)現(xiàn)，球軸承在高速重載條件下亞表面易產(chǎn)生微裂紋形成疲勞源，擴(kuò)展至表面形成剝落，進(jìn)而慢慢出現(xiàn)大面積表面剝落[9-12]。高速球軸承接觸部位會(huì)產(chǎn)生極高的接觸應(yīng)力和摩擦熱，循環(huán)往復(fù)的接觸應(yīng)力是球軸承疲勞失效的主要原因[13-15]。

對(duì)于軸承的接觸疲勞失效，國(guó)內(nèi)外學(xué)者從接觸表面剝落，微觀結(jié)構(gòu)演化，裂紋萌生及擴(kuò)展等方面進(jìn)行了大量研究：文獻(xiàn)[16]通過試驗(yàn)研究了接觸載荷對(duì)滾動(dòng)接觸損傷行為的影響，發(fā)現(xiàn)球與溝道接觸表面易萌生裂紋并向亞表面進(jìn)行擴(kuò)展；文獻(xiàn)[17]結(jié)合有限元法與斷裂力學(xué)理論對(duì)鋼軌滾動(dòng)接觸疲勞裂紋進(jìn)行應(yīng)力分析，并利用應(yīng)力強(qiáng)度因子描述裂縫尖端的應(yīng)力場(chǎng)；文獻(xiàn)[18]在循環(huán)載荷作用下建立模擬損傷的模型，得到循環(huán)載荷與損傷變量的關(guān)系，利用曲線擬合參數(shù)驗(yàn)證試驗(yàn)結(jié)果；文獻(xiàn)[19]采用擴(kuò)展有限元法與內(nèi)聚力單元相結(jié)合的方法研究了載荷加載方向?qū)?fù)合型裂紋擴(kuò)展的影響；文獻(xiàn)[20-22]基于概念數(shù)據(jù)模型的方法模擬鋼軌滾動(dòng)接觸中材料亞表層裂紋的萌生和擴(kuò)展。以上皆是在滾動(dòng)速度較小的工況下進(jìn)行的數(shù)值模擬和試驗(yàn)分析，大多將剪應(yīng)力視為疲勞壽命預(yù)測(cè)中起主導(dǎo)作用的失效應(yīng)力，結(jié)合數(shù)值計(jì)算對(duì)亞表面裂紋萌生和擴(kuò)展進(jìn)行模擬，并未研究在高速工況下摩擦生熱、材料夾雜硬度對(duì)數(shù)值分析結(jié)果的影響。

鑒于此，本文提出一種在高速工況下新的數(shù)值分析方法，設(shè)置可以熱力耦合的內(nèi)聚力單元，考慮摩擦熱和不同硬度的夾雜，并觀察其對(duì)整個(gè)微觀結(jié)構(gòu)的影響。首先，應(yīng)用循環(huán)迭代法通過Fortran語言對(duì)軸承徑向、軸向載荷進(jìn)行編譯，在Visual Studio上運(yùn)行模擬載荷并加載到有限元仿真軟件中，得到聯(lián)合載荷與應(yīng)力場(chǎng)的關(guān)系；其次，在球軸承與內(nèi)外圈溝道的接觸區(qū)域亞表面建立0厚度的內(nèi)聚區(qū)模型，設(shè)置具有溫度自由度的內(nèi)聚力單元；最后，編寫Python隨機(jī)插入圓形夾雜并引入摩擦生熱的溫度場(chǎng)進(jìn)行熱力耦合分析。通過對(duì)高速工況下球軸承接觸亞表面裂紋萌生的微觀研究，基于斷裂力學(xué)和損傷力學(xué)數(shù)值分析，將軸承微觀溫升數(shù)值結(jié)果與宏觀試驗(yàn)結(jié)果進(jìn)行比較，為高速球軸承接觸區(qū)域微觀機(jī)理分析提供模型支撐。

1 理論模型的構(gòu)建

深溝球軸承球與內(nèi)外圈接觸，其赫茲接觸狀態(tài)決定了接觸區(qū)域各點(diǎn)處的應(yīng)力、應(yīng)變狀態(tài)，且接觸應(yīng)力的分布呈高度的局域性。球與溝道接觸區(qū)域微觀模型如圖1所示，模擬了球與溝道接觸處的表面下方即亞表面處的接觸狀態(tài)，圖中：Fr為徑向載荷，F(xiàn)a為軸向載荷，P為載荷應(yīng)力。建立赫茲接觸狀態(tài)下的應(yīng)力應(yīng)變與內(nèi)聚力單元中的拉力分離之間的本構(gòu)關(guān)系，為使模型可以忽略黏結(jié)處的分層，選用基于拉力-分離(σ-δ)的0厚度內(nèi)聚力單元，通過自定義Python腳本實(shí)現(xiàn)在任意實(shí)體單元間對(duì)模型全局插入0厚度的內(nèi)聚力單元，如圖2所示。

圖1 球與溝道接觸區(qū)域微觀模型圖Fig.1 Micro model diagram of contact area between ball and raceway

圖2 實(shí)體單元間插入內(nèi)聚力單元Fig.2 Inserting cohesive elements between solid elements

用拉力-分離(σ-δ)方程代替工程中應(yīng)力-應(yīng)變(σ-ε)方程，內(nèi)聚力單元應(yīng)力傳遞關(guān)系如圖3所示，α為裂紋張開角度。

圖3 內(nèi)聚力單元應(yīng)力傳遞關(guān)系Fig.3 Stress transfer relation of cohesive elements

基于損傷力學(xué)與斷裂力學(xué)公式[23]，應(yīng)力與分離距離的關(guān)系為

σi=Kiδi，

式中：σi為表面應(yīng)力；Ki為材料剛度；δi為分離距離。

材料剛度與損傷變量的關(guān)系為

損傷初始準(zhǔn)則為

損傷起始后的分離距離為

在ABAQUS中，彈性剛度的退化可以通過損傷變量來表示。損傷變量(退化剛度SDEG)為

2 數(shù)值模擬

2.1 模型的建立

本文使用有限元分析軟件ABAQUS模擬高速球軸承接觸微觀區(qū)域，創(chuàng)建長(zhǎng)0.2 mm、寬0.1 mm的二維平面表示球與溝道接觸區(qū)域，在模型幾何中心靠近表面處進(jìn)行分區(qū)細(xì)化，設(shè)定分區(qū)長(zhǎng)為0.14 mm，寬為0.07 mm。網(wǎng)格類型為雙線性四節(jié)點(diǎn)的四邊形單元(CPE4R)，為提高計(jì)算收斂性，采用減縮積分和沙漏控制。中間分區(qū)部位網(wǎng)格單元設(shè)置為0.5 μm，其他部分網(wǎng)格大小由內(nèi)向外從0.5 μm到1 μm遞增，網(wǎng)格總數(shù)為18 207，節(jié)點(diǎn)數(shù)為36 911。球與溝道接觸區(qū)域的網(wǎng)格劃分如圖4所示。

圖4 球與溝道接觸區(qū)域的網(wǎng)格劃分Fig.4 Meshing of contact area between ball and raceway

模型材料參數(shù)見表1，內(nèi)聚力單元的材料參數(shù)見表2。

表1 模型材料參數(shù)Tab.1 Model material parameters

表2 內(nèi)聚力單元材料參數(shù)Tab.2 Material parameters of cohesive elements

2.2 聯(lián)合載荷下的接觸應(yīng)力

為模擬高速球軸承在接觸微觀區(qū)域的受力情況，首先應(yīng)用循環(huán)迭代法通過Fortran語言對(duì)軸承徑向、軸向載荷進(jìn)行編程，接著在Visual Studio中進(jìn)行模擬計(jì)算，最后加載到有限元仿真軟件建立的接觸區(qū)域模型中，即可得到聯(lián)合載荷作用下球軸承微觀接觸區(qū)的應(yīng)力場(chǎng)分布。

仿真中的載荷為壓力，壓力施加的起點(diǎn)為-0.02 mm，每一分析步壓力作用區(qū)域移動(dòng)的距離為0.006 mm，即球軸承運(yùn)轉(zhuǎn)速度為10 000 r/min，徑向壓力和軸向壓力的初始值設(shè)為P0=2 GPa。為探究徑向壓力與軸向壓力對(duì)亞表面處最大接觸應(yīng)力的影響，軸向壓力為初始值時(shí)，徑向壓力分別設(shè)置為1，2，3 GPa；徑向壓力為初始值時(shí)，軸向壓力分別設(shè)置為1，2，3 GPa。徑向、軸向壓力分別與亞表面最大接觸應(yīng)力之間的關(guān)系如圖5所示：隨著軸承轉(zhuǎn)數(shù)的累加，接觸應(yīng)力均略微上升；增大徑向壓力，接觸應(yīng)力也隨之較大幅度增大，增大軸向壓力，接觸應(yīng)力先增大后減??；兩者對(duì)比反映出徑向壓力對(duì)接觸應(yīng)力的影響占主導(dǎo)地位。

圖5 聯(lián)合載荷與接觸應(yīng)力的關(guān)系Fig.5 Relationship between combined load and contact stress

徑向壓力作用下接觸區(qū)域的應(yīng)力場(chǎng)分布如圖6所示：在接觸區(qū)域無損傷的情況下，最大接觸應(yīng)力在接觸表面下方即亞表面。這與赫茲接觸理論恰好吻合，為后續(xù)裂紋分析提供依據(jù)。在聯(lián)合載荷作用下，當(dāng)亞表面的內(nèi)部應(yīng)力大于材料的損傷極限應(yīng)力時(shí)便會(huì)產(chǎn)生疲勞損傷累積。

圖6 徑向壓力作用下接觸區(qū)域的應(yīng)力場(chǎng)分布Fig.6 Stress field distribution in contact area under radial load

2.3 夾雜與溫升對(duì)亞表面裂紋的影響

編寫Python腳本在微觀接觸模型亞表面內(nèi)聚區(qū)隨機(jī)插入半徑為r=3 μm、材料為硫化錳和氧化鋁的圓形夾雜，如圖7所示。硫化錳的彈性模量為138 GPa，泊松比為0.3；氧化鋁的彈性模量為375 GPa，泊松比為0.23。在模型上施加載荷進(jìn)行應(yīng)力分析，并開發(fā)可以進(jìn)行熱力耦合的內(nèi)聚力單元。在有限元模型文件里定義接觸摩擦生熱、對(duì)流換熱、熱傳導(dǎo)和輻射換熱，在接觸屬性中設(shè)置滾動(dòng)接觸的摩擦生熱及初始溫度，改變模型文件單元類型，設(shè)置可實(shí)現(xiàn)內(nèi)聚力單元傳熱的熱單元。

圖7 隨機(jī)插入的圓形夾雜Fig.7 Randomly inserted circular inclusions

當(dāng)損傷變量達(dá)到1時(shí)，不同夾雜的應(yīng)力云圖及亞表面處最大應(yīng)力位置變化如圖8所示，氧化鋁夾雜邊界剛度退化情況如圖9所示。

(a) 硫化錳夾雜

圖9 氧化鋁夾雜邊界剛度退化Fig.9 Boundary stiffness degradation of alumina inclusions

由圖8可知，引入夾雜后亞表面處的最大接觸應(yīng)力的位置與未引入夾雜時(shí)最大接觸應(yīng)力的位置相比更靠近接觸表面；對(duì)于硫化錳夾雜，亞表面最大接觸應(yīng)力位置距離接觸表面7.5 μm，向上偏移4.8 μm；對(duì)于氧化鋁夾雜，亞表面最大接觸應(yīng)力位置距離接觸表面7.7 μm，向上偏移4.6 μm；隨著深度的增加，應(yīng)力減小，夾雜材料對(duì)應(yīng)力值和應(yīng)力分布有著顯著影響。由圖9可知，損傷變量已達(dá)到1，根據(jù)斷裂力學(xué)此時(shí)裂紋開始出現(xiàn)；放大觀察夾雜邊界的剛度退化可以看出離亞表面越遠(yuǎn)的區(qū)域，裂紋萌生的概率越小；離亞表面越近的區(qū)域，裂紋萌生的概率越大。通過對(duì)不同材料夾雜的模擬，可看出接觸區(qū)域裂紋的萌生主要由亞表面最大應(yīng)力的位置及其向外衰減的趨勢(shì)所決定。將夾雜材料設(shè)為較軟的氧化鋁后，接觸區(qū)域的最大應(yīng)力數(shù)值有所下降，與材料較硬的夾雜相比，較軟的夾雜邊界不易萌生裂紋。

其他條件不變，引入溫度場(chǎng)進(jìn)行熱力耦合分析，摩擦生熱后的應(yīng)力云圖和熱流量云圖如圖10所示。摩擦生熱條件下，最大接觸應(yīng)力位置與接觸表面幾乎重合，并且夾雜邊界出現(xiàn)了一定程度的應(yīng)力集中。結(jié)合圖9與圖10可以看出，熱流與夾雜交匯處的邊界易形成裂紋且裂紋大多在夾雜的左右兩端。

(a) 熱應(yīng)力云圖

引入摩擦前后損傷變量隨軸承轉(zhuǎn)數(shù)的變化關(guān)系如圖11所示，在相同的軸承轉(zhuǎn)數(shù)下，考慮摩擦生熱影響的曲線會(huì)更快達(dá)到損傷變量1(損傷閾值)，即達(dá)到裂紋開裂狀態(tài)，因此摩擦熱對(duì)球軸承裂紋的萌生起促進(jìn)作用。

在聯(lián)合載荷作用下，接觸亞表面裂紋和夾雜附近應(yīng)力分布如圖12所示，接觸亞表面的裂紋和夾雜附近出現(xiàn)“蝴蝶”形應(yīng)力分布，這與已有研究[24]發(fā)現(xiàn)的起源于次表面夾雜物的蝶狀組織一致，說明了模型的正確性。

圖11 引入摩擦前后損傷變量隨軸承轉(zhuǎn)數(shù)的變化Fig.11 Variation of damage variable with bearing rotation before and after introducing friction

圖12 聯(lián)合載荷作用下接觸亞表面裂紋和夾雜附近應(yīng)力分布Fig.12 Stress distribution near contact subsurface cracks and inclusions under combined load

3 試驗(yàn)驗(yàn)證

為驗(yàn)證模型的有效性，進(jìn)行溫度試驗(yàn)，使用的T3-8N/2雙工位球軸承疲勞壽命試驗(yàn)機(jī)如圖13所示，主要由試驗(yàn)機(jī)主體、動(dòng)力系統(tǒng)、加載系統(tǒng)、潤(rùn)滑系統(tǒng)、電氣控制系統(tǒng)、計(jì)算機(jī)監(jiān)控系統(tǒng)等組成。采用608深溝球軸承，通過試驗(yàn)機(jī)的加載螺母加載徑向載荷和軸向載荷，循環(huán)105次。試驗(yàn)機(jī)采用計(jì)算機(jī)自動(dòng)控制，可自動(dòng)采樣與存儲(chǔ)數(shù)據(jù)。

圖13 高速球軸承疲勞壽命試驗(yàn)機(jī)Fig.13 High-speed ball bearing fatigue life tester

轉(zhuǎn)速為10 000 r/min時(shí)，軸向載荷為25 N，徑向載荷分別為50，100，150，200 N，軸承試驗(yàn)溫度與數(shù)值仿真模型溫度的變化如圖14a所示。轉(zhuǎn)速為10 000 r/min時(shí)，徑向載荷為25 N，軸向載荷分別為50，100，150，200 N，軸承試驗(yàn)溫度與數(shù)值仿真模型溫度的變化如圖14b所示。

圖14 軸承試驗(yàn)溫度與數(shù)值仿真模型溫度隨載荷的變化Fig.14 Variation of bearing test temperature and numerical simulation model temperature with load

由圖14可知，數(shù)值模型溫度略大于試驗(yàn)溫度。這是因?yàn)樵囼?yàn)檢測(cè)到的是軸承外圈溫度，無法檢測(cè)到內(nèi)圈溝道處溫度，故而造成溫度偏小。但就總體趨勢(shì)而言，試驗(yàn)數(shù)據(jù)與數(shù)值模型數(shù)據(jù)具有較高的一致性，反映出所建高速球軸承接觸微觀區(qū)域數(shù)值仿真模型的合理性。

4 結(jié)論

本文以高速工況下球軸承接觸微觀區(qū)域?yàn)檠芯繉?duì)象，將赫茲接觸與基于損傷力學(xué)的內(nèi)聚力單元本構(gòu)關(guān)系作為理論基礎(chǔ)，建立高速球軸承微觀接觸區(qū)域模型，分析球軸承接觸微觀區(qū)域亞表面應(yīng)力場(chǎng)的分布、夾雜與溫升對(duì)裂紋的影響，主要結(jié)論如下：

1)高速球軸承接觸微觀區(qū)域在聯(lián)合載荷作用下，徑向載荷對(duì)接觸應(yīng)力的影響占主導(dǎo)地位，接觸區(qū)域裂紋的萌生主要是由亞表面最大接觸應(yīng)力的位置以及其向外衰減的趨勢(shì)所決定。

2)夾雜使亞表面處最大接觸應(yīng)力的位置向上移動(dòng)，隨著深度的增加，應(yīng)力減小，剛度退化衰減；距離亞表面越遠(yuǎn)的區(qū)域裂紋萌生的概率較小，距離亞表面越近的區(qū)域裂紋萌生的概率較大；與材料較硬的夾雜相比，較軟的夾雜邊界不易萌生裂紋。

3)在摩擦生熱的條件下，夾雜邊界出現(xiàn)一定程度的應(yīng)力集中，熱流與夾雜交匯處的邊界上易形成裂紋且裂紋大多在夾雜的左右兩端，摩擦熱對(duì)球軸承裂紋的萌生起促進(jìn)作用。