考慮材料微觀結構的滾動接觸疲勞壽命研究

黃彥彥，占煜村，楊建平，駱孝武，祝晉旋

(1.成都大學 機械工程學院,四川 成都 610106;2.東方汽輪機有限公司,四川 德陽 618000;3.四川大學 空天科學與工程學院,四川 成都 610065)

0 引 言

作為旋轉機械的關鍵零部件，滾動軸承廣泛應用于人造衛星、空間機械臂等航天裝備中，其服役壽命對航天裝備整體服役性能具有重要影響.疲勞失效是滾動軸承主要失效形式[1]，對于安裝正確、潤滑良好的滾動軸承，其損壞形式一般為疲勞失效.根據Lundberg-Palmgren滾動接觸疲勞壽命理論[2]，接觸應力是滾動接觸疲勞壽命預測的重要參數.選取合適的材料模型，獲取更為準確的應力分布，對于疲勞壽命分析結果具有顯著影響.張若凡等[3]從疲勞強度、疲勞裂紋萌生與小裂紋擴展的微結構敏感性3個方面，對近年的研究成果進行了闡述，分析了金屬材料在微結構性質、形態及分布等特征下的超長壽命疲勞行為與失效機理.因此，開展材料表征方法對滾動接觸疲勞壽命預測影響的研究具有重要意義.

滾動軸承接觸疲勞失效研究一直以來都是國內外研究者關注的熱點問題.謝俊杰等[4]建立了考慮裂紋擴展的滾動接觸有限元分析模型，對裂紋的擴展過程進行了模擬分析.蔡森等[5]基于Hertz接觸理論，結合接觸疲勞壽命曲線，對不同載荷和速度下的高鐵軸承疲勞壽命進行了求解分析.金燕等[6]基于L-P疲勞壽命理論和Hertz接觸理論，開展了溫度對滾動軸承疲勞壽命影響規律的研究，得到了溫度對軸承壽命具有顯著影響的結論.孫玉鳳等[7]結合ANSYS和Fe-Safe軟件確定了滾動軸承最易發生疲勞失效的部位，并探究了殘余應力深度對疲勞壽命的影響規律.然而，上述研究均基于材料的各向同性假設，而忽略了金屬由許多細小的晶粒組成，在微觀上具有各向異性力學行為這一基本事實.近年來，相關學者們從模擬軸承材料真實微觀組織結構入手，開展滾動軸承疲勞壽命研究，得到了許多有益的成果.Raje等[8]利用Voronoi多邊形對軸承鋼多晶結構進行表征，分析了滾動軸承疲勞壽命具有統計學效應的微觀機理.Paulson等[9]基于Voronoi多邊形表征晶粒結構，結合彈流潤滑(EHL)和有限元仿真，對EHL狀態下的滾動軸承壽命進行了分析研究.Slack等[10]建立了基于Voronoi多邊形的顯式有限元模型，利用連續損傷理論對滾動軸承裂紋擴展行為進行了模擬仿真，并對軸承的壽命進行了預測.上述研究雖然引入了材料晶粒結構對滾動軸承接觸疲勞壽命進行分析，但均基于晶體彈性假設，未對晶粒塑性的影響進行考慮.

本研究從微觀角度出發，以Voronoi多邊形表征馬氏體軸承鋼晶粒結構，采用ABAQUS有限元軟件分別建立基于各向同性材料假設和考慮晶體學屬性的滾動接觸疲勞壽命分析模型，對兩種材料表征方法下滾動軸承相對疲勞壽命進行預測分析，探究晶體學屬性的滾動接觸疲勞壽命分析模型的特性.在此基礎上，開展晶粒尺寸分布對接觸疲勞壽命分布的影響研究，評估晶體學屬性滾動接觸模型的疲勞壽命預測能力.

1 建模方法

1.1 晶體塑性本構

(1)

(2)

(3)

式中，hαβ為硬化模量，/MPa；qαβ為描述潛硬化行為的矩陣；h0為初始硬化模量,/MPa；τs為飽和臨界分解剪切應力,/MPa；n為硬化系數.

1.2 晶體塑性有限元建模

1.2.1 赫茲接觸建模

滾子與滾道的接觸符合赫茲接觸方式.為降低模型規模，減少模型的計算時間，本研究將滾動軸承模型簡化為赫茲力加載的線接觸模型.本研究所建立有限元模型晶粒總數為5 000個，且平均晶粒大小為10 μm[8]，如圖1(a)所示.模型中各組成晶粒其取向各不相同，晶粒取向通過Python隨機數函數生成，共5 000種不同的晶粒取向，所建立有限元模型網格分布如圖1(b)所示.模型采用平面應變單元CPE4，共222 111個單元，223 112個節點.模型邊界條件根據軸承實際運行狀態進行設置，即對AB、BC、CD 3條邊自由度進行完全限制；AD邊為自由接觸表面，不限制其自由度.軸承鋼主要成分為馬氏體，本模型假設軸承鋼完全由馬氏體組成，不考慮殘余奧氏體等組成成分的影響.馬氏體力學性能參數如表1所示.滾道與滾動體的接觸載荷P如圖1(a)所示，可由式(4)計算得到，

表1 馬氏體力學性能參數

(a)有限元網格分布

(4)

式中，Pmax為最大接觸應力，為0.5 GPa；b為接觸半寬，設為25 μm.

軸承鋼晶粒尺寸發生變化，其強度和韌性將隨之改變.總體而言，對于晶粒尺寸大于100 nm的金屬，常溫下晶粒尺寸越小，其強度越高、疲勞性能越好[13].本研究以平均晶粒尺寸為10 μm的馬氏體軸承鋼為研究對象，引入晶體塑性開展滾動軸承疲勞壽命研究，并與基于各向同性材料假設的模型分析結果進行對比分析.彈性模量和屈服強度為各向同性材料模型建模時的必要參數，為準確確定上述材料參數，建立了平均晶粒尺寸為10 μm的馬氏體薄板拉伸試驗有限元分析模型.所建立有限元模型網格分布如圖2所示，模型共有85 330個單元，104 904個節點，單元類型為C3D8.采用一端固定，一端以1.125 μm/s的速度移動的方式進行拉伸試驗.

圖2 馬氏體多晶拉伸試驗模型

有限元分析得到的加載端力—位移曲線，經過一定的換算可以獲得相應的應力—應變曲線，如圖3所示.通過對應力—應變曲線彈性變形階段斜率的計算，得到馬氏體薄板彈性模量約為362 GPa，計算值與文獻[14]實驗測得值一致，說明了模型的有效性.由圖3可知，通過模型得到的應力—應變曲線沒有明顯的屈服點.因此，本研究采用工程上常用的σ0.2對材料屈服強度值進行估算，屈服強度約為1 350 MPa.

圖3 拉伸試驗應力—應變曲線

2 結果與討論

2.1 應力分布

圖4所示為Pmax為0.5 GPa時，基于各向同性材料假設和考慮軸承鋼晶體學屬性模型的下表面等效(Eguivalent)應力分布情況.由圖4可知，兩種情況的應力分布趨勢基本一致，最大von Mises應力均出現在下表面某處，而不是位于接觸表面.同時，基于各向同性材料假設的模型其應力過渡十分平滑，而對于考慮晶體學屬性的模型，其應力分布呈現出一定的離散性.由于組成材料各晶粒取向各不相同，當受力方向相同時，每個晶粒將表現出不同的應力狀態，最終導致整體應力表現出一定的離散狀態.此外，從圖4中還可以看出，考慮晶體學屬性的模型其最大von Mises應力值要高于基于各向同性材料假設的模型，這主要是由晶界處所產生的應力集中造成的.

(a)各向同性材料

圖5給出了von Mises應力沿赫茲接觸中心線AB、CD的變化情況.由圖5可知，對于基于各向同性材料假設的模型和考慮晶體學屬性的模型，兩者的von Mises應力值均呈現出迅速增大到最大值，然后逐漸減小的規律，且二者應力最大值發生位置基本重合.此外，AB、CD線上相應位置處的整體應力值并無明顯差別，說明考慮晶體學屬性的模型，整體上同樣滿足赫茲接觸應力分布.對于基于各向同性材料假設的模型，其von Mises應力變化十分平滑，而對于考慮晶體學屬性的模型，由于晶粒取向不同的影響，應力變化呈現出劇烈振動的現象.

圖5 赫茲接觸中心線von Mises應力

2.2 壽命分布

滾動軸承壽命存在多種預測方法，如Lundberg等[2]，Ioannides等[15]，Zhou[16]和Zaretsky[17]等提出的方法.其中，Lundberg等壽命預測模型應用最為廣泛，且其有效性已被大量實驗所驗證.因此，本研究選取Lundberg方法對滾動軸承接觸壽命進行評估.Lundberg滾動接觸疲勞壽命預測模型如下，

(5)

式中，S為材料生存概率；N為疲勞壽命循環數，/次；e為Weibull斜率；c為應力指數，τ為臨界剪切應力，/MPa;z為最大應力發生位置深度，/mm;V為體積應力，/MPa;A、c、h為由實驗確定的材料常數.對于軸承鋼，c和h分別取為10.33和2.33.由上式可知，生存概率S與臨界應力值成反比，與出現最大應力的深度成正比.因此，式(5)可簡化[8]為,

N∝zr/τq

(6)

式中，∝表示兩個變易為正比例關系.為方便起見，假設指數q和r與c和h相同，即r為2.33，q為10.33.

建立9個考慮晶體學屬性的赫茲接觸模型，對其應力分布進行提取，并利用式(6)對接觸疲勞壽命進行預測.為便于與基于各向同性材料假設的模型進行對比，本研究以相對壽命表征滾動軸承疲勞壽命.圖6為基于各向同性材料假設模型與考慮晶體學屬性模型預測的滾動軸承相對疲勞壽命.由圖6可知，對于考慮晶體學屬性的模型，其預測的軸承平均壽命不到基于各向同性材料假設模型的三分之一.由于基于材料各向同性假設的模型未考慮晶粒間的應力集中效應，所以其預測的疲勞壽命相較考慮晶體學屬性的模型偏大.從圖6中還可以看出，9個考慮晶體學屬性模型預測的相對疲勞壽命并不完全一致，這是由于材料具有不同的微觀結構，應力分布不均勻造成的.相較于基于各向同性材料假設的模型，考慮晶體學屬性的模型具有能夠分析壽命統計性分布的優點.

圖6 滾動接觸相對疲勞壽命

2.3 累積塑性剪切應變

對于各向同性材料模型，當外力不足以使材料產生大于屈服強度的應力時，塑性變形便不會發生.本研究的馬氏體屈服強度為1 350 MPa.顯然，采用Pmax=0.5 GPa不足以使馬氏體軸承鋼產生塑性變形.圖7(a)為Pmax=2.0 GPa時考慮晶體學屬性模型的累積剪切應變狀態.由圖可知，當Pmax為2.0 GPa時，材料內部會產生一定的塑性應變，最大累積剪切應變值為1.030×10-3(應變為形變量與原尺寸比，無量綱).當Pmax為0.5 GPa時考慮晶體學屬性模型的累積剪切應變狀態如圖7(b)所示.與基于各向同性材料假設的模型不同，即使外力大小未達到能使材料產生塑性變形的程度，材料內部仍有極小(相對于2.0 GPa的情況)的累積剪切應變產生.

(a)Pmax=2.0 GPa

綜上所述，相對于基于各向同性材料假設的模型，考慮晶體學屬性的模型能夠更好地對材料塑性變形行為進行分析.

2.4 晶粒尺寸分布對接觸疲勞的影響

晶粒的尺寸分布具有隨機性，但總體上符合正態分布規律.本研究將討論晶粒尺寸分別在N(0.01 mm,1×10-8mm2)、N(0.01 mm,2.5×10-7mm2)、N(0.01 mm,1×10-6mm2)、N(0.01 mm,2.5×10-5mm2)和N(0.01 mm,1×10-4mm2)5種正態分布情況下馬氏體鋼的疲勞壽命分布情況，以進一步說明考慮晶體學屬性的滾動接觸疲勞模型在分析疲勞壽命離散性上的優勢.

將上述5種晶粒尺寸分布下得到的21個臨界應力及位置深度帶入式(6)中，得到相應的相對壽命，并將其擬合得到Weibull壽命曲線，如圖8所示.由圖8可知，晶粒尺寸分布的方差越小，對應的Weibull斜率越大，馬氏體鋼的壽命分布范圍越小，說明晶粒的均勻分布程度與疲勞壽命密切相關.晶粒分布越不均勻，晶界處產生較大應力集中的可能性越大，在很大程度上將降低材料的疲勞壽命.

圖8 不同晶粒尺寸分布下馬氏體鋼的接觸疲勞壽命分布

3 結 論

1)考慮晶體學屬性的模型計算出的von Mises應力要大于基于各向同性材料假設的模型的應力，且von Mises應力分布的離散性要高于基于各向同性材料假設的模型應力分布的離散性.

2)考慮晶體學屬性的模型其預測的接觸疲勞壽命不到基于各向同性材料假設模型疲勞壽命的三分之一，且能夠對滾動軸承壽命分布的統計性進行模擬表征.

3)在外力引起的應力小于材料屈服強度時，考慮晶體學屬性的模型仍能產生一定的塑性變形.

4)晶粒尺寸分布離散性越大，威布爾斜率越小，即接觸疲勞壽命具有更大的分散性.