滾動軸承疲勞失效故障的數字孿生虛擬實體建模

顧偉，張文遠，王恒

(1.啟東開放大學，江蘇啟東 226200；2.南通大學機械工程學院，江蘇南通 226019)

0 前言

滾動軸承具有轉速高、結構緊湊、壽命長等優點，在工業生產中得到廣泛的應用。滾動軸承在使用的過程中產生的各種故障如果不能被及時發現和有效診斷，將會對工業生產造成重大安全隱患。如何有效診斷滾動軸承疲勞損傷程度和故障類型一直是故障診斷領域研究的重點問題[1]。

振動信號監測法是滾動軸承故障診斷的重要方法，目前常用的信號處理與故障特征提取的方法有短時傅里葉變換(STFT)、經驗模態分解(EMD)、經驗小波變換(EWT)等[2]。但是，現有的對滾動軸承振動信號進行時域頻域分析的故障診斷方法大多都存在著會受噪聲的干擾、不能適應復雜多變的工作環境、故障樣本數據難以獲取等不足。

滾動軸承自問世以來，人們對于其工作狀態和性能退化過程的機制研究也逐漸深入。如L-P模型[3]以及I-H模型[4]，分析了軸承套圈次表面裂紋的產生與發展過程，并由此發展出了滾動軸承的疲勞壽命模型，后來CHIU和TALLIAN發現滾動軸承還有另一種即表面裂紋發展而來的磨損失效，并提出了Chiu-Tallian工程模型[5]，為滾動軸承的疲勞定壽提供了新的依據；FUKATA等[6]提出了基于Hertz接觸理論的兩自由度滾動軸承動力學模型， SAWALHI和RANDALL[7]在此基礎上引入了單位諧振器，最終形成了五自由度滾動軸承動力學模型。但現有研究主要針對軸承套圈局部缺陷以及滾動體故障等軸承嚴重故障類型的建模，缺乏對滾動軸承早期故障的模型研究，尤其是針對軸承早期故障微弱信號難以有效提取故障信息。

近年來，數字孿生(Digital Twin)概念逐漸被越來越多的學者關注并使用。數字孿生以數字化的方式建立物理實體的多維、多時空尺度、多學科、多物理量的動態虛擬模型來仿真和刻畫物理實體在真實環境中的屬性、行為、規則，最早被應用在軍工及航空航天領域[8-9]。目前，數字孿生的應用逐漸向各個領域擴展。任濤等人[10]開展了基于數字孿生的機載光電探測系統性能退化建模研究，基于動態貝葉斯網絡(DBN)對系統性能退化建模；龐宇[11]開展了基于數字孿生的設備多部件成組預防性維護方法研究，建立了基于非線性Weiner退化過程的部件退化模型；方圓等人[12]開展了基于數字孿生的設備加工過程監測技術研究；葉倫寬[13]基于數字孿生進行了旋轉設備故障診斷技術研究，建立了旋轉設備的數字模型。以上研究成果表明：利用數字孿生建立研究對象的多維虛擬模型，并結合物理實體在真實環境中的行為開展研究的理念，可以廣泛應用于設備維護、異常檢測等各個領域。通過建立滾動軸承故障的數字孿生虛擬實體，與滾動軸承物理實體數據信息融合生成孿生數據，為滾動軸承故障診斷提供了一種新的思路。

疲勞失效是滾動軸承的常見失效形式，本文作者對滾動軸承內外圈局部缺陷模型進行了改進，改進的軸承內外圈局部缺陷模型同時考慮滾動體運行路徑和滾動體-滾道沖擊力的影響，并且提出了一種針對次表面裂紋早期故障的建模方法，彌補了滾動軸承早期微弱故障振動機制研究的不足。

1 滾動軸承疲勞失效故障數字孿生體

GRIEVES教授提出的數字孿生模型包括3個部分，即物理實體、虛擬實體以及二者之間的連接[14]。在不斷的實際應用中，這一框架得到了發展與補充，目前國內學者普遍認可的數字孿生框架是陶飛教授提出的一種數字孿生五維模型[8]，如式(1)所示：

MDT=(PE,VE,Ss，DD,CN)

(1)

式中：PE表示物理實體;VE表示虛擬實體;Ss表示服務，即模型所實現的功能；DD表示孿生數據；CN表示各組成部分間的連接。根據式(1)所示的數字孿生框架，數字孿生驅動的滾動軸承故障診斷流程如圖1所示。通過建立軸承故障的虛擬實體模型(VE)，并與滾動軸承物理實體(PE)相關參數信息融合形成孿生數據(DD)，利用融合后的孿生數據判斷軸承是否發生故障，并實現故障類型的診斷。

圖1 數字孿生驅動的滾動軸承故障診斷Fig.1 Rolling bearing fault diagnosis driven by digital twin

虛擬實體(VE)建模是數字孿生的重要組成部分，本文作者在五自由度滾動軸承振動模型的基礎上，利用改進的內外圈局部缺陷模型和滾道次表面裂紋機制模型建立滾動軸承疲勞失效故障的全壽命數字孿生虛擬實體，可以用來模擬疲勞失效引起的軸承故障從早期、中期到嚴重階段的發展歷程。通過文中提出的滾動軸承疲勞失效故障數字孿生虛擬實體，與滾動軸承物理實體數據信息融合生成軸承全壽命周期各個階段的孿生數據，為基于數字孿生的滾動軸承故障診斷提供了新的思路。

2 五自由度非線性滾動軸承局部故障機制模型

2.1 滾動軸承五自由度振動模型

文中所采用的是滾動軸承五自由度振動模型，主要研究內外圈的水平與豎直方向振動，將滾動體視為剛體，如圖2所示。該模型包括內外圈水平和豎直共4個自由度，以及單位諧振器1個自由度[7]。單位諧振器用來模擬受激軸承以及其他部件的高頻固有振動。

圖2 滾動軸承五自由度振動模型Fig.2 Vibration model of rolling bearing with five DOF

模型假設滾道與滾動體的接觸滿足Hertz接觸理論，各組成部分采用集中質量法，忽略各組件的轉動慣量與軸承滾道波紋度的影響。圖2中：ks、kp、kr分別為內圈、外圈及軸承座、單位諧振器的剛度；cs、cp、cr分別為內圈、外圈及軸承座、單位諧振器的阻尼；ms、mp、mr分別為內圈、外圈及軸承座、單位諧振器的質量；xs、xp、ys、yp、yb分別為內外圈水平與豎直兩個自由度的位移以及單位諧振器數值方向的位移。

非線性接觸力fx與fy是關于內圈轉動速度ωs和時間t的函數。軸承受純徑向載荷時，會在軸承滾道產生載荷區與非載荷區，載荷區中的滾動體會產生彈性變形使套圈之間產生相對位移，并產生柔性振動。

第j個滾動體的總變形量為

δj=(xs-xp)cosφj+(ys-yp)sinφj-cj=1,2,…,nb

(2)

(3)

(4)

式中：φj表示第j個滾動體的角位置；φ0表示滾動體初始角位置；ωc表示保持架角速度；ωs表示內圈角速度；Db表示滾動體直徑；Dp表示軸承節圓直徑；c表示游隙；nb表示滾動體個數。

根據赫茲接觸理論，第j個球與滾道之間的接觸力為

(5)

式中：kb為滾動體與內外圈之間總載荷-變形系數，單位為N/m1.5；δj為內外圈接觸變形。

當滾動體經過載荷區時，將式(2)代入式(5)可以得到單個滾動體與滾道之間的接觸力，累加得到套圈所受到的總的非線性接觸力在x與y方向上的分量fx與fy。

2.2 改進的內外圈局部缺陷故障模型

軸承故障可以分為套圈故障、滾動體故障和保持架故障，其中內外套圈的疲勞剝落失效和磨損失效是滾動軸承最常見失效形式[3]。如圖3所示，內外圈滾道產生單點局部缺陷時，當滾動體經過缺陷區，產生的變形與經過無缺陷區域時會發生變化，其總變形量為

δj=(xs-xp)cosφj+(ys-yp)sinφj-c-βjcd

(6)

式中：cd為局部缺陷深度；βj為開關函數，當滾動體經過缺陷區時，βj值為1，否則βj值為0；c表示游隙；φj表示第j個滾動體的角位置。當軸承外圈存在局部缺陷時，設局部缺陷跨越的角度為Δφd，缺陷的初始角度為φd，則當φd<φj<φd+Δφd時βj的值為1，否則βj值為0。軸承內圈存在局部缺陷時，缺陷位置會隨著內圈轉動而轉動，此時，缺陷的初始角度為

φd=ωst+φd0

(7)

式中：ωs表示內圈轉動角速度；t表示時間；φd0表示時間t=0時的缺陷初始角位置。

圖3 外圈局部缺陷模型示意Fig.3 Local defect model of outer ring

現有的局部缺陷模型僅根據缺陷產生的尺寸不同和Hertz接觸理論指出，滾動體經過缺陷時產生的接觸變形會服從不同的分段函數[15]，但未考慮滾動體經過缺陷區邊緣時撞擊滾道產生的脈沖沖擊對振動響應的影響[16]。本文作者綜合考慮接觸力與沖擊力，提出一種改進的內外圈局部缺陷模型。當滾動體經過缺陷時，滾動體到達的缺陷深度隨著球的旋轉而逐漸增加。缺陷深度達到最大值時，與滾道發生撞擊，隨后深度逐漸減小，在撞擊發生的時刻，滾動體與滾道會產生接觸力與沖擊力。所以，滾動體位移與轉動角度的函數為

(8)

滾動體所受合力與轉動角度的函數為

(9)

式中：φ為滾動體進入缺陷區后轉過的角度，取值范圍為[0, Δφd]；fc為因撞擊產生的沖擊力；kb為滾動體與內外圈之間總載荷-變形系數；δj為內外圈接觸變形；cd為局部缺陷深度。

2.3 滾道次表面裂紋早期故障模型

根據L-P模型[3]，滾動軸承套圈滾道的疲勞裂紋首先出現在最大交變剪應力處，也就是靠近滾道表面的次表面處，然后向表面擴展。這種疲勞裂紋的產生，是軸承服役周期的重要階段，這種次表面疲勞裂紋不同于第2.2節的局部缺陷模型，不會影響滾動體的運行路徑，其接觸力仍然滿足Hertz接觸理論，但由于次表面裂紋的存在，總載荷-變形系數kb和接觸變形會產生變化。文獻[17]研究了滾動軸承的剛度特性，軸承剛度與載荷、接觸變形的關系近似為

Q=kδ

(10)

式中：Q表示載荷；k表示接觸剛度；δ表示接觸變形。從式(10)可以看出，當滾動軸承發生次表面裂紋早期故障，有滾動體經過故障區時其接觸剛度會發生變化。

外圈存在滾道次表面裂紋故障時：

(11)

內圈存在滾道次表面裂紋故障時：

G2(ks)=

(12)

式中：kp為外圈接觸剛度；ks為內圈接觸剛度；λ是與接觸變形量有關的故障特征參數，可以根據帶有次表面裂紋的滾動軸承受載時，內外圈產生的變形計算得出。

3 滾動軸承全壽命數字孿生虛擬實體模型構建

根據L-P模型[3]，內外套圈上產生的疲勞裂紋引發的剝落是影響軸承使用壽命的顯著因素。全壽命周期模型將滾動軸承故障的產生及發展分為3個階段，即在軸承潤滑良好、載荷正常的情況下，軸承長期服役會產生疲勞失效，由于受滾道與滾動體的Hertz接觸力影響，疲勞裂紋首先出現在靠近軸承內外圈滾道表面一定深度的次表面處，次表面裂紋產生為早期故障階段；次表面裂紋會向表面擴展，逐漸形成表面裂紋和疲勞剝落失效形式，表面裂紋和疲勞剝落產生為中期故障階段；表面裂紋與疲勞剝落兩種故障的出現，會導致軸承周期性振動幅度增大，軸承摩擦增大，故障區不斷發展最終導致內外圈滾道發生大范圍的嚴重磨損失效，使軸承嚴重失效甚至卡死停止工作，此時為嚴重故障階段。疲勞失效下的滾動軸承全壽命數字孿生虛擬實體模型如圖4所示。

圖4 滾動軸承全壽命數字孿生虛擬實體模型

3.1 滾動軸承早期故障虛擬實體模型

次表面裂紋的出現被稱為滾動軸承早期故障階段，如第2.3節所述，次表面裂紋的出現會導致軸承剛度特性的變化。通過建立滾動軸承的有限元模型可以計算出早期故障軸承的接觸變形。以6220軸承為例，在內外圈滾道受最大剪應力處去除一部分軸承體積單元來模擬次表面裂紋的產生，有限元幾何模型如圖5所示。對滾道施加Hertz接觸力，經過靜力學仿真得出滾道表面的接觸變形如圖6所示。

圖5 次表面裂紋幾何模型

圖6 接觸變形仿真結果

不同尺寸的次表面裂紋所產生的接觸變形如圖7所示，軸承受徑向載荷50 kN，當不存在次表面裂紋時，可以看出隨著次表面裂紋的擴大，接觸變形也會隨之增大。根據公式(10)，當作用在套圈的總載荷不變時，可以通過不同尺寸的次表面裂紋所產生的接觸變形與不帶有裂紋時產生的接觸變形的比值得出第j個滾動體經過早期故障區時的故障特征參數λ的值，將λ代入公式(11)、(12)可得滾動軸承次表面裂紋早期故障模型。

圖7 次表面裂紋長度與接觸變形的關系

3.2 滾動軸承中期故障虛擬實體模型

當軸承進入中期故障階段。根據Chiu-Tallian工程模型[11]，表面裂紋可能是由次表面裂紋發展而來，也可能由于次表面最大剪應力有時會小于材料的疲勞極限，疲勞裂紋會直接出現在滾道表面。內外圈次表面裂紋與表面裂紋一旦出現，就會不斷發展，形成滾道表面的小范圍剝落。不同于次表面裂紋，表面裂紋與疲勞剝落會影響滾動體經過故障區時的運行路徑，所以文中將表面裂紋與疲勞剝落稱為中期故障階段。

當表面裂紋產生時，可用第2.2節介紹的內外圈局部缺陷模型，此時故障跨越角度遠小于故障最大深度，滾動體不會接觸到缺陷的最大深度，其運行軌跡如圖8(a)所示。滾動體經過缺陷產生的位移與轉動角度的曲線如圖8(b)所示。

圖8 內外圈表面裂紋早期故障模型

表面裂紋會不斷擴大導致內外圈滾道的疲勞剝落，根據第2.2節介紹的內外圈局部缺陷模型，疲勞剝落產生的缺陷區跨越角度與缺陷最大深度尺寸相近，滾動體經過缺陷區時，滾動體到達缺陷深度不斷增大，隨后與滾道缺陷產生一次沖擊，隨后深度逐漸減小。滾動體的運行軌跡如圖9(a)所示，滾動體經過缺陷產生的位移與轉動角度的曲線如圖9(b)所示。

圖9 內外圈疲勞剝落中期故障模型

3.3 滾動軸承嚴重故障虛擬實體模型

內外圈表面裂紋會發展成為疲勞磨損，疲勞磨損故障相比疲勞剝落故障具有更大的缺陷區跨越角度，同時，故障區對滾動體運行路徑的影響也更加復雜。此時，滾動體經過缺陷區時，到達缺陷區的深度不斷增大并到達最大深度，同時產生一次沖擊，滾動體在缺陷區運行一段時間后，產生第二次沖擊，同時到達缺陷區深度逐漸減小。滾動體的運行軌跡如圖10(a)所示，滾動體經過缺陷產生的位移與轉動角度的曲線如圖10(b)所示。

圖10 內外圈疲勞磨損中期故障模型

4 滾動軸承故障數字孿生虛擬實體模型結果分析

將圖4中的各種故障形式按文中介紹的方法建立動力學模型，然后分別做數值仿真，形成軸承全壽命周期的振動信號虛擬數據，由此建立滾動軸承全壽命周期數字孿生虛擬實體模型(VE)。通過數值仿真方法，在時間域內對系統響應的時間歷程進行離散，把運動微分方程分為各離散時刻的方程，并進行逐步數值積分，求解出軸承系統一系列離散時刻上的瞬態響應值。文中采用中心差分振動仿真方法。中心差分法將系統的運動微分方程化成對時間的差分格式，利用逐步積分求出一系列離散時刻的響應值。

系統在t時刻的運動方程為

(13)

選取6220深溝球軸承作為研究對象，對所提出的滾動軸承故障模型進行振動仿真分析。對軸承施加50 kN的徑向載荷，軸承外圈固定，內圈轉速為1 500 r/min。6220深溝球軸承的幾何參數如表1所示，運動微分方程相關參數如表2所示。

表1 6220深溝球軸承相關參數Tab.1 6220 deep groove ball bearing related parameters

表2 軸承動力學微分方程相關參數[13]

針對早期故障階段和中期故障階段中的外圈帶有局部缺陷、內圈帶有局部缺陷、外圈帶有次表面裂紋3種情況進行研究。其中，內、外圈局部缺陷的初始角度φd為0°，缺陷跨越角度Δφd為3°，外圈局部缺陷的最大深度為1.5 mm，內圈局部缺陷的最大深度為0.5 mm，此時滾動體運行路徑如圖10所示，會發生兩次沖擊。次表面裂紋缺陷跨越角度Δφd為5°。根據L-P模型，次表面裂紋首先在最大交變剪應力處產生， 6220軸承受載荷時的剪應力平均深度最大交變剪應力出現在深度0.75 mm處，這表明了裂紋向表面擴展時，深度最大可達到1.5 mm。此次仿真設定次表面裂紋長度為1 mm，缺陷區跨越角度為5°，根據圖7，故障特征參數λ為0.96。

4.1 內外圈局部缺陷模型數據結果分析

對外圈帶有局部缺陷、內圈帶有局部缺陷2種軸承故障模型進行數值仿真，得到振動仿真信號如圖11、圖12所示。觀察圖11(a)與圖12(a)可以得到：內圈或外圈出現局部缺陷時的時域信號是一系列脈沖衰減運動，這是由于滾動體與滾道產生了周期性的沖擊振動。外圈的沖擊振動周期為0.009 2 s，內圈沖擊振動周期為0.012 7 s。對比圖13與圖12(a)可知：改進內、外圈局部缺陷模型相比傳統內、外圈局部缺陷模型可以更有效地模擬缺陷-滾動體沖擊響應，與文獻[17]所述的沖擊響應特征一致性較好。觀察圖11(b)可以得到：外圈局部缺陷頻譜信號在112.69 Hz以及它的二倍頻、三倍頻附近具有較高的譜線，且幅值隨著頻率的增加而衰減；觀察圖12(b)可以得到：內圈局部缺陷頻譜信號的主頻率為75.49 Hz，與內圈周期0.0127 s的倒數相近，且幅值隨著頻率的增加而衰減。當滾動軸承外圈上有單個損壞點時，其包絡譜圖是一系列離散的譜線，這些線幅值隨著頻率的增加而逐漸減小，以外圈缺陷的特征頻率為間隔；而內圈的頻譜以內圈缺陷為特征的譜線在每個階次的多個頻率處振幅逐漸減小。與仿真的結果相符，從理論上驗證了局部缺陷模型和振動仿真結果的準確性。

圖11 外圈局部缺陷虛擬模型振動信號

圖12 內圈局部缺陷虛擬實體模型振動信號

圖13 不考慮內外圈局部缺陷仿真信號

4.2 次表面裂紋模型數據結果分析

圖14(a)為外圈滾道出現次表面裂紋早期故障的振動仿真信號時域波形，圖14(b)為它的頻譜圖。觀察圖14(a)可以得到:早期故障的時域波形不同于局部缺陷時域信號，不是周期性的脈沖衰減運動，當滾動體經過故障區時，軸承的振動響應幅值會有所增大，但最大幅值遠小于局部缺陷振動響應。觀察圖14(b)可以得到：早期故障頻譜信號的特征頻率與外圈缺陷頻譜信號相似，幅值隨著頻率的增加而衰減。

圖14 外圈滾道次表面裂紋虛擬實體振動信號

5 結論

建立了滾動軸承全壽命周期數字孿生虛擬實體模型，為實現基于數字孿生的滾動軸承故障診斷提供了模型基礎?；跀底謱\生的滾動軸承故障診斷可以克服傳統的數據驅動和基于機制模型的故障診斷方法受噪聲干擾大、不能適應復雜多變的工作環境、故障樣本難以獲取等缺點。文中主要研究成果如下：

(1) 在滾動軸承五自由度非線性振動模型的基礎上，考慮滾動體-滾道沖擊力的影響，改進了軸承局部缺陷模型，并基于載荷-接觸變形關系建立了次表面裂紋引起的滾道次表面裂紋早期故障模型。

(2) 根據疲勞失效理論，對滾動軸承疲勞失效的產生與發展歷程建模，建立了滾動軸承全壽命周期數字孿生虛擬實體模型。全壽命模型將軸承生命周期分為早期故障、中期故障和嚴重故障階段，可以用來表征軸承故障的產生與發展到嚴重失效的過程。

(3) 利用有限元法確定了次表面裂紋故障模型特征參數。對次表面裂紋故障模型和改進內、外圈局部缺陷模型進行數值仿真，形成虛擬模型振動信號。對虛擬模型數據結果進行分析，結果表明：兩種故障模型信號的故障頻域特征與理論相近，且改進內、外圈局部缺陷模型可以更有效地模擬滾道缺陷-滾動體的沖擊響應，驗證了文中提出的故障模型的正確性。