基于二分法-PSO的碰摩轉(zhuǎn)子有限元模型修正

劉 濤，劉 娟

（1.呼和浩特鐵路局 包頭車輛段，內(nèi)蒙古 包頭 150200；2.中國石油大學(xué)（北京）機械與儲運工程學(xué)院，北京 102249）

旋轉(zhuǎn)機械在我們的生活中起著非常重要的作用，伴隨著現(xiàn)代工業(yè)的發(fā)展，旋轉(zhuǎn)設(shè)備向著高精度，高速化和智能化的方向發(fā)展，從而導(dǎo)致轉(zhuǎn)子和定子間的間隙越來越小，結(jié)構(gòu)越來越緊密，轉(zhuǎn)子碰摩故障發(fā)生得越來越頻繁[1]。

轉(zhuǎn)子與定子碰摩是一類非常復(fù)雜的動力學(xué)問題。采用剛?cè)狁詈系霓D(zhuǎn)子有限元模型，研究一個或少數(shù)幾個參數(shù)對碰摩響應(yīng)行為的影響，為常用的研究路徑[2]。如何建立一個精度較高的含有碰摩故障的轉(zhuǎn)子模型，對轉(zhuǎn)子碰摩故障的研究，是很重要的。以模態(tài)試驗為依據(jù)，進行模型修正，可以有效地解決模型精度較低的問題。從1968年Fox和Kapoor首次提出特征向量靈敏度的計算方法開始，基于靈敏度的模型修正方法得到廣泛應(yīng)用[3]。徐張明，高天明和Amar Bouazzouni，Mohammad等都對靈敏度模型修正方法進行了研究，且取得了一定的成果[4–6]。現(xiàn)在應(yīng)用較為廣泛的方法為基于有限元模型的設(shè)計參數(shù)型修正方法，其優(yōu)點是物理意義清晰，工程實用性更強。其中，基于靈敏度分析的設(shè)計參數(shù)型模型修正方法應(yīng)用最廣泛[7]。隨著工程需求的日趨提高，人們對模型修正技術(shù)的計算精度和計算效率也有了更高的要求。

為了提高剛?cè)狁詈限D(zhuǎn)子模型柔性轉(zhuǎn)軸的精度，以及模型修正技術(shù)的計算精度與計算效率。本文基于模態(tài)試驗數(shù)據(jù)，將二分法–PSO作為模型修正的優(yōu)化算法，實現(xiàn)有限元模型修正，從而提高了有限元模型修正的速度與精度，且避免了在SPO迭代尋優(yōu)過程中陷入局部最優(yōu)解的問題。基于模型修正理論，利用轉(zhuǎn)子系統(tǒng)模態(tài)試驗獲取的轉(zhuǎn)軸固有頻率，以及WorkBench建立的有限元模型，進行模型修正。研究結(jié)果表明，修正后的有限元模型精度有了大幅度的提高，可以較為準確地反應(yīng)試驗?zāi)Ｐ偷慕Y(jié)構(gòu)特征。最后利用修正后的模型進行轉(zhuǎn)子碰摩故障動力學(xué)特性的研究，考慮轉(zhuǎn)速變化對碰摩故障的影響，通過分析結(jié)果可以知道修正后的轉(zhuǎn)子模型可以較為準確地反映試驗?zāi)Ｐ偷倪\行結(jié)果。

1 模型修正方法理論

1.1 靈敏度模型修正

以靈敏度為基礎(chǔ)的迭代修正法與有限元模型的離散本質(zhì)一致，只需要修正系統(tǒng)中的某幾個單元或節(jié)點的參數(shù)，就可以使有限元模型迅速收斂于與試驗數(shù)據(jù)良好相關(guān)的修正模型，因此備受青睞[8–10]。目標函數(shù)是最優(yōu)的設(shè)計性能，它是有關(guān)設(shè)計變量的函數(shù)

其中p代表設(shè)計參數(shù)向量，{fe}，{fp(p)}分別代表結(jié)構(gòu)動態(tài)特性的實驗值和分析值；r(p)為殘差項；VLB、VUB分別代表結(jié)構(gòu)設(shè)計參數(shù)變化的上下限；Wf代表結(jié)構(gòu)各個特征量之間的加權(quán)矩陣。

靈敏度的計算公式為

則修正公式可以寫成

其中S?為靈敏度矩陣S的廣義逆，因此只要得到系統(tǒng)的靈敏度矩陣，確定了特征量的變化量就可以通過優(yōu)化迭代法確定參數(shù)的變化量。

1.2 粒子群優(yōu)化算法

PSO是一種隨機搜索全局優(yōu)化新方法，相較于其他算法，具有實現(xiàn)容易、精度高、收斂快等特點[11]。

對于有N個粒子群Xi（xi1，xi2，…，xid）的速度，位置更新公式為

其中Vi(vi1，vi2，…，vid)為始速度，Pbesti為粒子Xi找到的最優(yōu)點，Gbest為所有粒子找到的最優(yōu)點，C1，C2為加速度因子，r1，r2是均值，用以增加粒子的多樣性，ω為慣性權(quán)。

1.3 二分法–PSO

PSO算法在迭代過程中隨機性較強，要達到高精度需要較長時間。為此將二分法引進迭代中，由盲目搜索變成計劃搜索。

因此，按照圖1中的流程，完成對轉(zhuǎn)子系統(tǒng)的模型修正。

圖1 基于二分法–SPO的模型修正總體思路

具體算法：

(1)初始化X，V，Gbest，Pbest；

(2)更新X，V；

(3)給定的搜索區(qū)間：

計算函數(shù)Fj=f(X1min，X2min，…，Xjmax，…Xnmax)，

其中j=1，2…n；選取最小值Fi=min(F1，F(xiàn)2，…Fn)，即Fi=f(X1min，X2min，…，Ximax)；

(4)按式(5)或者式(6)更新參數(shù)X；

(5)判斷：設(shè)定誤差e，若Fi＜e，進行步驟6，若若Fi＞e，返回步驟3；

(6)更新Pbest，Gbest；

(7) 判 斷 如 果 iter＞iterMax，或 者 Gbest＞GbestMax，則退出計算，否則返回步驟2。

2 分析模態(tài)

2.1 試驗?zāi)B(tài)分析

為了獲得試驗數(shù)據(jù)，通過對自由-自由狀態(tài)下的轉(zhuǎn)子結(jié)構(gòu)進行激振，測量激振力和響應(yīng)數(shù)據(jù)，進行數(shù)據(jù)采集和信號分析，從而獲得模態(tài)數(shù)據(jù)[12]。試驗時用充氣的輪胎墊住轉(zhuǎn)子結(jié)構(gòu)以模擬自由-自由狀態(tài)。如圖3所示。

圖2 模態(tài)試驗設(shè)計圖

圖3 轉(zhuǎn)子有限元模型

選擇轉(zhuǎn)軸的8個節(jié)點作為測點，采用單點激勵輸入、多點輸出方式，對轉(zhuǎn)子系統(tǒng)的每一個測點施加一定的激勵，用加速度傳感器采集激勵力信號，采集節(jié)點上的振動響應(yīng)信號。其中設(shè)備包括：轉(zhuǎn)子振動實驗臺，傳感器，力錘，等。根據(jù)采集到的激勵信號和響應(yīng)信號，識別出轉(zhuǎn)軸前4階固有頻率，如表1所示。

表1 轉(zhuǎn)軸前4階固有頻率

2.2 計算模態(tài)分析

在SolidWorks軟件中建立轉(zhuǎn)子系統(tǒng)的三維實體模型，模型構(gòu)建過程中忽略該結(jié)構(gòu)件中某些局部結(jié)構(gòu)，如圓角、倒角等細節(jié)以減少有限元計算時間，提高有限單元的協(xié)調(diào)性。將適當(dāng)簡化后的三維模型導(dǎo)入ANSYS Workbench進行有限元分析。由于模型中的軸承為滑動軸承。因此，通過彈性阻尼連接方式的Bushing連接進行模擬，由于阻尼較小忽略不計，而軸承剛度通過K值來計算。

邊界條件為：在轉(zhuǎn)軸的兩端以及中間部分添加圓柱約束，同時，在轉(zhuǎn)軸中安裝轉(zhuǎn)子的部分添加應(yīng)力。網(wǎng)格劃分：采用Workbench自動網(wǎng)格劃分技術(shù)，智能地控制網(wǎng)格大小，使得計算精度不會太低，又可以保證計算速度不會太低。最后進行受力分析，建立有限元模型，同時進行模態(tài)分析計算，得到相應(yīng)的模態(tài)參數(shù)。轉(zhuǎn)子系統(tǒng)有限元模型如圖3所示，圖4為轉(zhuǎn)軸前4階模態(tài)的振型圖。

3 模型修正結(jié)果

參數(shù)的選擇與試驗中靈敏度有很大的關(guān)系。首先考慮可能影響有限元模型的幾個特征參量：彈性模量E，密度M，質(zhì)量系數(shù)M、C，轉(zhuǎn)軸的截面半徑R，轉(zhuǎn)軸的長度L，轉(zhuǎn)軸的支承剛度K。而各個參數(shù)對模型特性的影響程度不同，鑒于參數(shù)靈敏度是選擇優(yōu)化參數(shù)的一個重要指標，表示各參數(shù)對模型特性的影響程度。因此，可以通過靈敏度的分析確定需要優(yōu)化的參量。

圖4 轉(zhuǎn)軸前4階模態(tài)特性

如圖5所示，各參數(shù)在各階次模態(tài)下的靈敏度是不同的。為了提高迭代優(yōu)化的效率，防止在迭代過程中發(fā)生病理現(xiàn)象，在參數(shù)選擇中，需選擇靈敏度較高的參量。質(zhì)量系數(shù)MC與長度L的靈敏度相對于其他參量較小可以忽略不計，轉(zhuǎn)軸的截面半徑R，彈性模量E，軸承剛度K靈敏度相對較大，可以選擇為待修正的設(shè)計參量。同時，密度M的靈敏度雖不是很大，但對計算結(jié)果影響較大，因此也作為優(yōu)化參量。

圖5 優(yōu)化參數(shù)及其靈敏度分析

在模型修正的過程中，各階次的固有頻率與試驗結(jié)果間的誤差是不斷變化的。如圖6所示。

這是分別采用PSO算法與二分法–PSO算法對目標函數(shù)進行優(yōu)化迭代的過程。PSO算法與二分法–PSO算法都是伴隨著迭代次數(shù)的增加，目標函數(shù)值不斷變小，直到最后趨于穩(wěn)定。然而PSO算法對目標函數(shù)的優(yōu)化迭代在20次以上誤差逐漸穩(wěn)定且趨于0，而二分法–PSO方法對目標函數(shù)的優(yōu)化迭代在10次左右，誤差就逐漸穩(wěn)定且趨于0。同時，二分法–PSO方法每次迭代計算所用時間比PSO算法每次迭代計算所用時間少。經(jīng)過簡單計算可以看出，改進后的方法迭代計算所用時間少于傳統(tǒng)方法所用時間。

由此可以證明采用二分法–PSO算法對目標函數(shù)的優(yōu)化迭代的速度更快。同時得到表2所示的修正前后的設(shè)計參數(shù)值。

圖6 目標函數(shù)收斂過程

從表3的結(jié)果可以看出，修正前固有頻率與試驗分析獲取的固有頻率間的差異相對較大，說明仿真前的模型在反映試驗?zāi)Ｐ偷慕Y(jié)構(gòu)特性時，誤差相對較大。而修正后的有限元模型的計算結(jié)果與試驗測試結(jié)果更相近，即采用PSO方法與基于二分法–PSO方法進行修正后的結(jié)果都比初始模型分析結(jié)果的準確度要高，而且基于二分法–PSO方法修正結(jié)果的誤差比PSO方法修正的誤差更小。說明了修正方法及設(shè)計參數(shù)優(yōu)化結(jié)果的正確性，同時可以表明基于二分法–PSO的修正方法可以得到更高的計算精度。

4 模型應(yīng)用

圖7 轉(zhuǎn)子碰摩試驗

以圖7所示的試驗臺為仿真研究的原型，進行轉(zhuǎn)子系統(tǒng)碰摩故障的研究，其中碰撞位置采用塑料塊與轉(zhuǎn)子圓盤摩擦模擬碰摩故障。

將修正后的有限元模型導(dǎo)入ADAMS中，構(gòu)建剛?cè)崤己系霓D(zhuǎn)子模型，然后進行轉(zhuǎn)子系統(tǒng)碰摩故障的動力學(xué)仿真研究。由于試驗轉(zhuǎn)子的1階臨界轉(zhuǎn)速為5 160 r/min，偏心量為m?r=2.5×10-4kg?m，碰摩間隙設(shè)置為0.15 mm。在ADAMS中，碰摩故障的研究主要利用IMPACT函數(shù)來模擬兩個構(gòu)件之間的碰撞力，IMPACT函數(shù)構(gòu)造的碰撞力由彈性力與阻尼力組成的，由于在試驗中，用塑料塊與轉(zhuǎn)子圓盤摩擦模擬碰摩故障，因此根據(jù)ADAMS推薦，碰摩剛度為1 150 N/m，阻比系數(shù)為0.68，滲透深度為0.1，碰撞指數(shù)為2。然后添加約束條件，即轉(zhuǎn)子與轉(zhuǎn)軸間的固定副，軸與底座間的旋轉(zhuǎn)副。添加一個恒定轉(zhuǎn)速，大小為9 000(度每秒)，即1 500 r/min來模擬電機的驅(qū)動，在轉(zhuǎn)子系統(tǒng)的另一邊添加一個負載轉(zhuǎn)矩0，為了降低突變帶來的誤差，因此使用STEP函數(shù)控制，使得速度穩(wěn)定增加到9 000(度每秒)。

圖8為修正前，修正后，實驗測量獲取的振動頻率圖。

從圖中可以看出，修正后模型的振動頻率與試驗測量的振動頻率都在26 Hz左右，而修正前的振動頻率在35 Hz附近，對比試驗結(jié)果，誤差較大，因此可以看出修正后的結(jié)果更為接近試驗結(jié)果。

圖10是轉(zhuǎn)速為2 760 r/min時的頻譜圖與軸心軌跡圖，通過圖9與圖10可以看出，在轉(zhuǎn)速低于臨界轉(zhuǎn)速時，存在大于基頻的高倍頻成分，軸心軌跡有內(nèi)凹現(xiàn)象，伴隨著轉(zhuǎn)速的增長，二倍頻增加的速度遠大于一倍頻增長的速度。

表2 修正前后的模型參數(shù)

表3 模型修正前后固有頻率與誤差的變化情況

圖8 碰摩轉(zhuǎn)子修正前后及實驗測量頻率圖

圖9 臨界轉(zhuǎn)速以內(nèi)的仿真結(jié)果

圖10 臨界轉(zhuǎn)速以內(nèi)的仿真結(jié)果

圖11為臨界轉(zhuǎn)速情況下的振動情況，此時轉(zhuǎn)速在5 160 r/min附近，出現(xiàn)小于一倍頻的成分，基頻幅值遠大于之前，同時，軸心軌跡為多環(huán)嵌套；圖12中，轉(zhuǎn)速在6 190 r/min附近，是臨界轉(zhuǎn)速的1.2倍，因此出現(xiàn)了1X/1.2的低頻成分，同時軸心軌跡為多環(huán)嵌套。

圖11 臨界轉(zhuǎn)速時仿真結(jié)果

圖12 臨界轉(zhuǎn)速以外仿真結(jié)果

5 結(jié)語

為了提高剛?cè)狁詈限D(zhuǎn)子模型柔性轉(zhuǎn)軸的精度，以及模型修正技術(shù)的計算精度與計算效率，本文利用模態(tài)試驗數(shù)據(jù)，應(yīng)用二分法–PSO算法實現(xiàn)有限元模型修正的優(yōu)化計算，從而提高了有限元模型修正的速度與精度，且避免了在SPO迭代尋優(yōu)過程中陷入局部最優(yōu)解的問題。對比PSO模型修正的結(jié)果與二分法–PSO模型修正的結(jié)果，發(fā)現(xiàn)基于二分法–PSO的模型修正方法在計算速度，計算精度方面相都有所提高。

將修正后的模型應(yīng)用在轉(zhuǎn)子碰摩故障的動力學(xué)特性研究中，考慮轉(zhuǎn)速變化對轉(zhuǎn)子碰摩動力學(xué)特性的影響。研究結(jié)果表明，當(dāng)轉(zhuǎn)速小于臨界轉(zhuǎn)速時，伴隨著轉(zhuǎn)速的增加二倍頻幅值的增長迅速；在臨界轉(zhuǎn)速時，振動幅值遠大于先前；當(dāng)轉(zhuǎn)速大于臨界轉(zhuǎn)速時，出現(xiàn)低于一倍頻的成分。仿真結(jié)果與試驗結(jié)果一致性較高，為實際工程中轉(zhuǎn)子系統(tǒng)的安全穩(wěn)定運行提供一定的實際參考。

[1]MING ZENG,YU YANG,JINDE ZHENG.Normalized complex Teager energy operator demodulation method and its application to fault diagnosis in a rubbing rotor system[J].Mechanical Systems and Signal Processing,2015(50-51):380-399.

[2]FERNANDO L PONTA,ALEJANDRO D OTERO,LUCAS I LAGO,Effects of rotor deformation in windturbine performance:The dynamic rotor deformation bladeelementmomentum model(DRD-BEM)[J].Renewable Energy,2016(92):157-170.

[3]FOX R L,KAPOOR M P.Rate of change of eigenvectors and eigenvalues[J].AIAA Journal,1968,12(6):2426-2429.

[4]徐張明，高天明，沈榮瀛，等.一種改進的利用頻響函數(shù)進行有限元模型修正的方法[J].振動與沖擊，2002，21(3)：43-45.

[5]ASMA F,BOUAZZOUNI A.Finite element model updating using variable separation[J].European Journal of Mechanics-A/Solids,2007,26(4):728-735.

[6]MOHAMMAD R D,JAVAD V A,SIROUS G,et al.Determination of nonlinear behavior of a ball joint system by model updating[J].Journal of Constructional Steel Research,2012,71:52-62.

[7]GUILLAUME ROSENZVEIG,FRAN?OIS LOUF,A FE model updating method for the simulation of the assembly process of large and lightweight aeronautical structures[J].Finite Elements in Analysis and Design,2016(111):56-63.

[8]BAIYONG GOU,WEIJIE ZHANG,QIUHAI LU.A successive selection method for finite element model updating[J].Mechanical Systems and Signal Processing,2016(70-71):320-333.

[9]NEHETE D V,MODAK S V,GUPTA K.Coupled vibroacoustic modelupdating using frequency response functions[J].Mechanical Systems and Signal Processing,2016(70-71):308-319.

[10]NYSSEN F,GOLINVAL J C.Identification of mistuning and model updating of an academic blisk based on geometry and vibration measurements[J].Mechanical Systems and Signal Processing,2016(68-69):252-264.

[11]ZHANG C,SUN J.An alternate two phases particle swarm optimization algorithm for flow shop scheduling problem[J].Expert Systems with Applications,2009,36(3):5162-5167.

[12]慶光蔚.基于復(fù)模態(tài)和頻響修正方法的磁軸承-轉(zhuǎn)子系統(tǒng)參數(shù)識別[D].南京：南京航空航天大學(xué)，2010.