劉旭南, 趙麗娟, 付東波, 張飛飛
(遼寧工程技術(shù)大學(xué) 機(jī)械工程學(xué)院,遼寧 阜新 123000)
采煤機(jī)是煤礦開采過程中最常見的綜采設(shè)備,其主要由截割部、牽引部、行走部和附屬裝置等部分組成,其中截割部消耗功率約占整機(jī)功率的90%左右,也是發(fā)生故障最高的部分[1-2],因此故障診斷模型的建立以及故障的有效診斷,對于采煤機(jī)維護(hù)、自動(dòng)化生產(chǎn)以及無人工作面的開發(fā)等方面顯得至關(guān)重要。
為此,專家學(xué)者也展開了相關(guān)研究。20世紀(jì)80年代,許景崑[3]以SIRUS-400型采煤機(jī)為研究對象,給出了系統(tǒng)的診斷方法和診斷標(biāo)準(zhǔn);20世紀(jì)90年代后,許多學(xué)者綜合了多種常見的分析方法,在采煤機(jī)故障現(xiàn)象采用正交小波變換法、模塊化神經(jīng)網(wǎng)絡(luò)和專家系統(tǒng)結(jié)合的方法開發(fā)了采煤機(jī)故障的智能專家系統(tǒng)[4-7];2006年楊高宏[8]引入斷齒故障模型,并結(jié)合動(dòng)力學(xué)仿真對齒輪嚙合力的頻率譜密度分析,確定設(shè)備故障情況下齒輪嚙合的頻率特征,進(jìn)而分析設(shè)備故障狀態(tài);2014年弓曉鳳等[9]用小波包變換對采煤機(jī)故障信號(hào)進(jìn)行分解重構(gòu)和能量計(jì)算,并以此提出一種減速器故障診斷方法。但是針對該方法的實(shí)驗(yàn)故障判斷不明確,不能精確判斷故障的發(fā)生部位;2015年劉曉樂等[10]通過研究不同齒輪斷齒形式下接觸力信息,分析斷齒故障形式對齒輪各階傳動(dòng)的影響,最終發(fā)現(xiàn)斷齒故障形式會(huì)使接觸力頻率邊頻帶能量變得密集;2017年郝志勇等[11]用銷軸傳感器對采煤機(jī)搖臂惰輪軸進(jìn)行工作載荷測量,提出一種相對準(zhǔn)確的惰輪軸載荷測量方法。本文采用MATLAB,PRO/E,ANSYS和ADAMS,建立包含斷齒、裂紋齒輪和缺陷軸承等常見故障的采煤機(jī)截割部截割煤巖的剛?cè)狁詈咸摂M樣機(jī)模型,通過仿真獲取各惰輪軸軸向和徑向受力數(shù)據(jù),通過小波包分解法將其轉(zhuǎn)化的各子帶能量值作為神經(jīng)網(wǎng)絡(luò)的輸入信號(hào),通過訓(xùn)練后,形成了虛擬樣機(jī)與神經(jīng)網(wǎng)絡(luò)相結(jié)合的小波包分解法采煤機(jī)截割煤巖的故障診斷模型。該模型可對復(fù)雜工況下服役的采煤機(jī)傳動(dòng)系統(tǒng)中齒輪裂紋、斷齒,軸承內(nèi)外圈、滾動(dòng)體、保持架缺陷等故障形式進(jìn)行預(yù)測和診斷,解決了以往僅能對單一故障零件的振動(dòng)信號(hào)進(jìn)行分析及只適用于分析截割普通煤巖工況的不足,對在線實(shí)時(shí)監(jiān)測具有一定的理論價(jià)值。
根據(jù)收集的常見故障案例[12],采用參數(shù)化建模方法建立傳動(dòng)系統(tǒng)的齒輪零件,根據(jù)裂紋擴(kuò)展軌跡[13]對齒輪1、齒輪2和行星輪引入裂紋特征和斷齒特征,以齒輪1為例,建立好的故障模型如圖1所示。
圖1 齒輪故障三維模型Fig.1 Gear failure three dimensional models
參照凱斯西儲(chǔ)大學(xué)2013年公布的軸承實(shí)驗(yàn)數(shù)據(jù)[14],建立寬度、深度分別為1.27 mm,0.7 mm的溝痕缺陷內(nèi)圈和外圈模型,建立直徑為0.18 ~0.043 mm、深度為0.28 mm的深坑點(diǎn)蝕滾動(dòng)體模型,相應(yīng)模型如圖2(a)~圖2(c)所示,而保持架缺陷為限制連續(xù)兩個(gè)滾動(dòng)體轉(zhuǎn)動(dòng),可通過Adams中對軸承保持架和滾動(dòng)體間的約束來實(shí)現(xiàn),其模型如圖2(d)所示。
柔性體可以更好模擬零部件工作時(shí)的接觸、及振動(dòng)等特征,因此將實(shí)驗(yàn)?zāi)P椭械臏y試軸、故障零件等關(guān)鍵零件做成柔性體形式。將PRO/E中建立好的三維模型直接導(dǎo)入到ANSYS中,按照設(shè)計(jì)要求定義材料屬性和實(shí)常數(shù),定義網(wǎng)格單元類型為solid 8node 185以模擬彈性接觸,并對齒輪、軸承的故障部位進(jìn)行局部細(xì)化;通過Rigid Region將零件工作中與其他零件的接觸面和該面幾何中心上的一點(diǎn)建立包含所有自由度的剛性區(qū)域;最后通過Export to ADAMS命令生成各零件對應(yīng)的MNF文件。
圖2 軸承故障三維模型Fig.2 Bearing failure three dimensional models
在PRO/E中建立采煤機(jī)截割部零件模型,并將其裝配成整體,通過MECH/pro導(dǎo)入ADAMS中,導(dǎo)入后的虛擬樣機(jī)模型如圖3所示。
圖3 虛擬樣機(jī)模型Fig.3 Virtual prototype model
在ADAMS中點(diǎn)選需要被替換的剛性零件,通過Rigid To Flex命令完成相關(guān)柔性零件的導(dǎo)入和替換,由于故障診斷系統(tǒng)以測試軸1(齒輪2軸)、測試軸2(齒輪2軸)、測試軸3(齒輪5軸)和測試軸4(齒輪8軸)的軸向和徑向受力信號(hào)為診斷參數(shù),故將其替換成對應(yīng)的柔性體文件,以便虛擬樣機(jī)仿真結(jié)束后提取數(shù)據(jù)。在替換后根據(jù)截割部實(shí)際情況添加約束、驅(qū)動(dòng)、接觸以及載荷,最終形成剛?cè)狁詈咸摂M樣機(jī)模型如圖4所示。其中比較重要的三個(gè)部分分別為:載荷和接觸的計(jì)算以及仿真步長的設(shè)定。
圖4 剛?cè)狁詈咸摂M樣機(jī)模型Fig.4 Rigid flexible coupling virtual prototype model
1.3.1 載荷的計(jì)算
采煤機(jī)載荷主要是由于截齒與煤巖相互作用產(chǎn)生的,如圖5所示,第j截齒在截割煤巖時(shí)會(huì)受到截割阻力Zj和牽引阻力Yj及側(cè)向力Xj[15],利用 MATLAB計(jì)算各個(gè)截齒的瞬態(tài)載荷,并將各個(gè)截齒不同的瞬態(tài)負(fù)載進(jìn)行轉(zhuǎn)化和合成,得到滾筒質(zhì)心處的載荷[16]。設(shè)置滾筒直徑φ800 mm,轉(zhuǎn)速83.48 r/min,牽引速度為4 m/min,煤巖堅(jiān)固性系數(shù)f=3,通過計(jì)算得到的載荷曲線如圖6所示。在ADAMS滾筒質(zhì)心點(diǎn)建立一個(gè)與圖5所示坐標(biāo)系一致的廣義力通過spline命令將其加載到該廣義力上以模擬滾筒載荷。
圖5 截齒受力示意圖Fig.5 Force diagram of the pick
圖6 載荷曲線Fig.6 Load curves
1.3.2 接觸的計(jì)算
接觸的添加最主要要準(zhǔn)確的計(jì)算接觸的剛度,對于無故障模型其剛度可用式(1)計(jì)算[17]
(1)
式中:kb1,kb2分別為主、從齒輪的彎曲剛度,N/mm,其計(jì)算如式(2)所示;kf1,kf2分別為主、從齒輪的輪體變形剛度,N/mm,其計(jì)算如式(3)所示;kh為赫茲接觸剛度,N/mm,其計(jì)算式(4)所示。
(2)
(3)
(4)
式中:E為彈性模量,MPa;W為齒輪齒寬,mm;α1為齒頂漸開線與齒形中線夾角的余角,(°);α2為齒形漸開線在基圓上所成弧夾角的一半,(°);α為壓力角,(°);L,M,P,Q為相關(guān)系數(shù)可通過式(5)獲得;uf為初始嚙合點(diǎn)曲率線和齒形中心線交點(diǎn)到齒根圓的距離,mm;sf為齒根圓上的齒形寬,mm;E1,E2為兩齒輪的彈性模量,MPa;R1,R2為接觸點(diǎn)距兩齒輪回轉(zhuǎn)中心的半徑,mm;v1,v2為兩齒輪的泊松比。
(5)
式中:hf=rf/r,rf為齒根圓半徑,r為齒輪輪轂內(nèi)徑,mm;X分別為L,M,P和Q,X為不同參數(shù)時(shí),A,B,C,D,E和F由表1確定。
表1 式(5)的系數(shù)
將齒輪1、齒輪2和行星輪相關(guān)參數(shù)其代入式(1)~式(5)便可獲得無故障狀態(tài)下齒輪模型的剛度:齒輪1與齒輪2嚙合剛度為964 463.85 N/mm;齒輪2之間的嚙合剛度為1 019 825.45 N/mm;行星齒輪與齒圈內(nèi)嚙合剛度為977 015.26 N/mm,行星齒輪與太陽輪外嚙合剛度714 241.82 N/mm。
對于故障區(qū)域以有限元法對其嚙合剛度進(jìn)行計(jì)算更加合理,通過對“1.2”節(jié)中故障齒輪進(jìn)行材料賦予、網(wǎng)格劃分添加約束和齒面受力F后,便可求解,在后處理中可以輸出齒面各節(jié)點(diǎn)的變形量δ,計(jì)算嚙合兩齒輪輪齒總變形量δ1和δ2,再根據(jù)嚙合總變形量Σδ計(jì)算輪齒剛度,總變形量Σδ,其值為δ1與δ2之和,剛度K和受力F之間的關(guān)系如式(6)所示。
(6)
式中:S為嚙合面積,mm。
通過計(jì)算得到在故障狀態(tài)下齒輪模型的剛度:齒輪1與齒輪2嚙合剛度為506 796.11 N/mm;齒輪2之間的嚙合剛度為581 302.05 N/mm;行星齒輪與齒圈內(nèi)嚙合剛度為428 544.6 N/mm,行星齒輪與太陽輪外嚙合剛度576 432.85 N/mm。
1.3.3 仿真步長的設(shè)定
截割部齒輪傳動(dòng)中最大振動(dòng)頻率應(yīng)為高速端齒輪嚙合頻率,其計(jì)算公式為
fc=aZ1n
(7)
式中:a為嚙合齒輪對數(shù);Z1為齒輪齒數(shù);n為轉(zhuǎn)速,r/s。
經(jīng)計(jì)算,高速端最大嚙合頻率fmax=974 Hz,根據(jù)采樣頻率定理,確定采樣頻率fs=2 000 Hz,故仿真步長為0.000 5 s。分別建立表2所示的11組仿真模型,并按相關(guān)參數(shù)設(shè)置完成虛擬仿真。
表2 仿真分組
在對11組模型進(jìn)行仿真后,對虛擬樣機(jī)試驗(yàn)所得的數(shù)據(jù)進(jìn)行提取,以測試軸1為例,從t=0.201 s開始,提取0.719 s(一滾筒轉(zhuǎn)動(dòng)周期)的數(shù)據(jù),由于仿真步長為0.000 5 s,故采樣點(diǎn)數(shù)為1 438,通過ADAMS后處理plot命令,獲取滾筒該轉(zhuǎn)動(dòng)周期正常和齒輪1斷齒工況下的測試軸1徑向力數(shù)據(jù),如圖7所示。
圖7 測試軸1徑向力Fig.7 Radial force of test shaft 1
其他仿真結(jié)果與圖7類似,但通過圖7無法辨識(shí)齒輪是否故障,因此可以采用數(shù)據(jù)處理手段將各組中的信號(hào)特征分離開,以便于故障診斷。
按照頻率差通過Coif小波包分解將各波峰分解到不同頻率子帶上,設(shè)置分解尺度為4,小波包分解各自帶頻率跨度為125 Hz,可以將各波峰分開。
對上節(jié)測試軸徑向力數(shù)據(jù)進(jìn)行小波包分解,可以得到16個(gè)子帶數(shù)據(jù),圖8僅顯示了其中的第1、第3、第5子帶。
圖8 小波包分解子帶Fig.8 Subbands of wavelet packet decomposition
采用相鄰兩點(diǎn)距離和作為該信號(hào)的能量值,通過式(8)獲取經(jīng)小波包分解后的各子帶信號(hào)能量計(jì)算信號(hào)能量,并將其進(jìn)行記錄,齒輪1斷齒情況下各測試軸的徑向和軸向小波包子帶能量值如表3所示。
(8)
Elman神經(jīng)網(wǎng)絡(luò)是一種動(dòng)態(tài)回歸型神經(jīng)網(wǎng)絡(luò),其網(wǎng)絡(luò)結(jié)構(gòu)如圖9所示。圖中矩形所代表的為承接層,每一個(gè)承接層節(jié)點(diǎn)都與一個(gè)隱層節(jié)點(diǎn)相連,存儲(chǔ)隱層回傳來的信息,并以權(quán)重w3傳給隱層,其值不隨時(shí)間變化,則有Elman神經(jīng)網(wǎng)絡(luò)數(shù)學(xué)模型。
表3 齒輪1斷齒下小波包子帶能量值
圖9 Elman神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.9 Thestructure of Elman neural network
(9)
式中:x為神經(jīng)網(wǎng)絡(luò)輸入值;y為神經(jīng)網(wǎng)絡(luò)輸出值;xy為隱層輸出量;xc為承接層輸出量;f(·)為隱層傳遞函數(shù);g(·)為輸出層傳遞函數(shù);b1,b2分別為隱層和輸出層的閾值;w1和w2分別為輸入層到隱層的權(quán)重和隱含層到輸出層的權(quán)重,隨迭代次數(shù)調(diào)整,w3為承接層到隱層權(quán)重,其值為定值不隨迭代次數(shù)變化。Elman神經(jīng)網(wǎng)絡(luò)評測值為實(shí)際輸出和目標(biāo)輸出的均方差
(10)
式中:C為均方差;n為迭代次數(shù);y(t)和yk(t)分別為網(wǎng)絡(luò)實(shí)際輸出和目標(biāo)輸出。
根據(jù)故障類型設(shè)置對應(yīng)各種故障類型的神經(jīng)網(wǎng)絡(luò)輸出向量,如表4所示。
表4 神經(jīng)網(wǎng)絡(luò)各工況對應(yīng)輸出值
以“2.2”節(jié)處理后得到的小波包子帶能量值數(shù)據(jù)為基礎(chǔ),將每次虛擬樣機(jī)仿真數(shù)據(jù)依照測試軸1、測試軸2、測試軸3和測試軸4的軸向、徑向16組子帶能量排列成128×1的列向量,作為神經(jīng)網(wǎng)絡(luò)訓(xùn)練樣本向量。采用MATLAB以不同隱層數(shù)、不同節(jié)點(diǎn)數(shù)建立不同結(jié)構(gòu)的Elman神經(jīng)網(wǎng)絡(luò),進(jìn)行多次實(shí)驗(yàn),確定精度最好的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)。最終選取如圖10所示的20-13-11三層Elman神經(jīng)網(wǎng)絡(luò),網(wǎng)絡(luò)訓(xùn)練迭代次為41次,訓(xùn)練精度可達(dá)7.06×10-7,該網(wǎng)絡(luò)對故障形式的分類收斂性較好,如圖11所示。
圖10 故障診斷神經(jīng)網(wǎng)絡(luò)Fig.10 Neural network for fault diagnosis
圖11 網(wǎng)絡(luò)均值誤差Fig.11 Mean squared error of the network
以MG2×70/325采煤機(jī)開采某礦區(qū)17層煤為例,按文獻(xiàn)[18]的方法在采煤機(jī)左截割部兩根齒輪2軸、齒輪5軸和齒輪8軸外壁上分別選取5個(gè)位置點(diǎn)安裝軸向和環(huán)向應(yīng)變片,將應(yīng)變片構(gòu)成橋式電路與外圍系統(tǒng)集成。圖12為采集該系統(tǒng)采集得到的軸向力、徑向力信號(hào)。
通過小波包分解后獲得的小波包子帶能量值如表5所示,其中測試軸1除軸向第4、第6、第8、第16子帶及徑向第9子帶與齒輪1斷齒下小波包子帶能量量級(jí)一致外,其余均均比齒輪1斷齒能量要大;測試軸2軸向單數(shù)子帶與齒輪1斷齒能量量級(jí)一致,而雙數(shù)子帶均比齒輪1斷齒能量要小1~2個(gè)量級(jí);測試軸2徑向3子帶、4子帶、7子帶、11子帶、12子帶、15子帶與齒輪1斷齒能量級(jí)要大一個(gè)量級(jí),其余子帶能量量級(jí)一致;測試軸3軸向和徑向子帶能量均要比齒輪1斷齒相應(yīng)能量要大1~2個(gè)量級(jí);而測試軸4與齒輪1斷齒能量基本類似。通過對比可知測試數(shù)據(jù)與齒輪1斷齒不屬于同一類故障,應(yīng)進(jìn)一步分析。而對于常見的11種工況可通過對比能量數(shù)據(jù)加以區(qū)別。
圖12 軸向力及徑向力信號(hào)Fig.12 Axial and radial force signals
子帶號(hào)測試軸1軸向徑向測試軸2軸向徑向測試軸3軸向徑向測試軸4軸向徑向1630.71 3 270 834.8 25.34 31 113.00 18.37 28 144.53 4 781.27 6 478 368.80 2213.78 471 576.05 23.12 17 140.10 17.49 26 266.42 683.49 771 252.25 3214.01 623 069.24 24.86 17 569.19 14.98 25 413.03 1 285.54 970 563.96 4233.70 606 451.04 21.98 18 144.39 15.57 26 741.69 723.39 919 867.00 5263.42 870 974.27 27.91 20 173.34 14.47 32 293.05 926.65 918 136.74 6215.54 555 511.70 21.59 16 692.61 16.29 26 937.99 722.00 970 816.86 7246.37 389 327.65 21.22 16 305.78 16.07 31 602.55 737.36 622 711.30 8213.88 417 619.88 22.67 21 166.12 17.63 29 546.55 637.48 671 720.85 9275.62 564 637.38 30.77 31 787.56 14.96 52 086.72 571.11 502 356.66 10204.45 329 288.43 27.05 15 460.14 15.54 30 195.99 356.30 362 063.25 11209.63 305 072.32 20.21 19 082.06 16.40 29 457.95 307.36 395 634.35 12206.75 294 839.74 25.86 18 003.54 15.55 25 112.54 317.18 376 480.49 13284.27 908 204.39 23.35 18 573.49 18.17 40 679.08 814.18 964 794.62 14204.89 438 150.03 24.36 16 001.10 14.97 33 002.88 466.27 676 800.73 15183.72 304 570.60 25.19 23 328.21 14.75 30 170.11 310.07 367 462.17 16244.47 432 273.73 24.25 19 070.6414.78 26 462.15 340.06 631 119.02
將表5數(shù)據(jù)整合成128×1的列向量輸入到建立好的Elman故障診斷神經(jīng)網(wǎng)絡(luò)中,得到的輸出向量為[0.050 5, 0.000 1, 0.000 2, 0, 0.159 2, 0.992 5, 0, 0, 0, 0.043 9, 0]T,與表4對比可知采煤機(jī)截割部此時(shí)為故障狀態(tài)下工作,故障形式為行星輪斷齒。這與現(xiàn)場實(shí)際情況相一致,如圖13所示,該采煤機(jī)截割截割部的減速器行星輪輪齒發(fā)生大面積斷裂,可見該方法可以有效診斷出截割部傳動(dòng)系統(tǒng)故障。
圖13 行星齒輪斷齒Fig.13 Planetary gear breakage
(1)綜合了多年來齒輪和軸承的故障特性研究和故障診斷方面的先進(jìn)方法,根據(jù)故障形式下齒輪和軸承的振動(dòng)特性,計(jì)算不同工況的采煤機(jī)滾筒載荷以及齒輪嚙合剛度,生成故障齒輪時(shí)頻嚙合剛度。
(2)建立剛?cè)狁詈咸摂M樣機(jī)模型,根據(jù)計(jì)算所得參數(shù)進(jìn)行虛擬樣機(jī)仿真,將仿真結(jié)果與理論值對比,證明仿真結(jié)果準(zhǔn)確可靠,并得到測試軸仿真結(jié)果數(shù)據(jù)。
(3)開發(fā)出一套以采煤機(jī)虛擬樣機(jī)實(shí)驗(yàn)為基礎(chǔ)提供設(shè)備運(yùn)行樣本數(shù)據(jù),以小波包變換對數(shù)據(jù)進(jìn)行特征提取,根據(jù)Elman神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)一種多參數(shù)采煤機(jī)搖臂故障診斷方法,使設(shè)備在投入使用之前就具有一套以虛擬樣機(jī)實(shí)驗(yàn)為樣本的適應(yīng)自身的故障診斷系統(tǒng)。
該研究對復(fù)雜工況下服役的采煤機(jī)故障檢測以及在線實(shí)時(shí)監(jiān)測具有一定的應(yīng)用價(jià)值。