周期時變轉(zhuǎn)子系統(tǒng)參數(shù)識別技術(shù)研究

王 志， 王建軍， 劉 玉

(1.北京航空航天大學 能源與動力工程學院，北京 100191； 2.沈陽航空航天大學 遼寧省航空推進系統(tǒng)先進測試技術(shù)重點實驗室，沈陽 110136)

隨著科技的進步，在航空航天、交通運輸、能源電力等領(lǐng)域中的工程結(jié)構(gòu)都在向復(fù)雜化、高速化和智能化的方向發(fā)展，這就導(dǎo)致許多結(jié)構(gòu)的參數(shù)是隨時間變化的，如火箭發(fā)射過程中的質(zhì)量、齒輪嚙合剛度等，這種時變特性常常使系統(tǒng)的動態(tài)特性發(fā)生重大變化，甚至會威脅到系統(tǒng)的穩(wěn)定性，因此時變系統(tǒng)的參數(shù)識別研究具有廣泛而重要的意義[1]。

目前，對于線性時不變結(jié)構(gòu)動力學正問題和反問題的研究已趨于成熟，但時變結(jié)構(gòu)的研究，因為其難度較大，還是結(jié)構(gòu)動力學學科的前沿問題[2]，特別是反問題(時變參數(shù)識別問題)的研究，被認為是近期力學領(lǐng)域的重點研究方向。因此國內(nèi)外一些學者對其進行了深入探索，并取得了豐碩的成果。史治宇、許鑫等利用系統(tǒng)的激勵和受迫響應(yīng)數(shù)據(jù)，提出了基于狀態(tài)空間和小波變換的時變系統(tǒng)參數(shù)識別方法，并以二自由度彈簧-質(zhì)量-阻尼模型為仿真算例，對系統(tǒng)時變參數(shù)進行了識別，算例結(jié)果驗證了其方法的正確性和有效性[3-6]。龐世偉等[7-10]提出了利用自由響應(yīng)數(shù)據(jù)的遞推子空間方法，并通過移動質(zhì)量塊-簡支梁模型的仿真試驗，驗證了該方法可以有效地辨識時變系統(tǒng)參數(shù)變化，且跟蹤精度高，計算量小，對噪聲不敏感。Shan等[11]利用連續(xù)小波變換(CWT)對時變系統(tǒng)進行了參數(shù)識別，并討論了小波尺度對識別結(jié)果的影響。Bao等[12]基于新改進的希爾伯特-黃變換理論(HHT)算法對參數(shù)時變的復(fù)合梁進行了有效識別，且魯棒性好。

以上都是以一般時變系統(tǒng)為研究對象，給出了仿真算例。在實際工程中存在著大量具有周期時變參數(shù)的振動系統(tǒng)(參激振動系統(tǒng))，如齒輪系統(tǒng)[13-14]、非對稱轉(zhuǎn)子(非圓截面軸、帶嵌線槽電機轉(zhuǎn)子等)或裂紋轉(zhuǎn)子系統(tǒng)等[15-16]，當前針對這種特定參激振動系統(tǒng)的參數(shù)識別研究還很少。

本文利用方波脈沖函數(shù)具有脫關(guān)性、正交性等優(yōu)良性質(zhì)，推導(dǎo)出參數(shù)識別遞推計算公式；并利用該公式分別對二自由度周期時變仿真系統(tǒng)和典型非對稱(方軸)轉(zhuǎn)子實驗系統(tǒng)進行了周期時變剛度識別。

1 參數(shù)識別理論及算法

結(jié)構(gòu)動力學問題，包括激勵、結(jié)構(gòu)本身和響應(yīng)三要素，即輸入、系統(tǒng)和輸出。參數(shù)識別問題，是指已知系統(tǒng)的輸入和輸出，求得描述系統(tǒng)特性的各種參數(shù)，屬于結(jié)構(gòu)動力學第一類逆問題。

一般線性時變n自由度系統(tǒng)由下述狀態(tài)方程描述

(1)

式中：x(t)∈Rn×1為狀態(tài)向量，A(t)∈Rn×n為時變系統(tǒng)矩陣，B(t)∈Rn×r為時變控制矩陣，u(t)∈Rr×1為輸入向量。參數(shù)識別即在已知x(t)、u(t)的情況下，識別出矩陣A(t)、B(t)中的各項參數(shù)。

本文基于系統(tǒng)在初始狀態(tài)下的自由響應(yīng)數(shù)據(jù)來識別A(t)，因此設(shè)u(t)=0。

1.1 函數(shù)的方波脈沖級數(shù)展開

一個在區(qū)間[0,T)內(nèi)絕對可積的函數(shù)f(t)可以在該時域內(nèi)展開成一組方波脈沖級數(shù)[17]，即

f(t)≈f1g1(t)+f2g2(t)+…+

(2)

式中：fi為第i項方波脈沖系數(shù)，gi(t)為第i項方波脈沖函數(shù)。其中g(shù)i(t)定義為

(3)

根據(jù)式(3)，經(jīng)推導(dǎo)可知，gi(t)具有脫關(guān)性與正交性，即

(4)

(5)

式(2)兩端同時乘以gi(t)，并在[0,T)上積分，由gi(t)的正交性質(zhì)可得

(6)

取其三次代數(shù)精度近似，則有

f(ih)]

(7)

其截斷誤差為

(8)

將式(1)中的x(t)、x0、A(t)分別展開為方波脈沖級數(shù)，且設(shè)N>n，則有

(9)

簡寫成

(10)

同理，可將x0、A(t)寫成

(11)

(12)

1.2 時變系統(tǒng)矩陣A(t)識別

當u(t)=0時，將式(1)中的狀態(tài)方程等號兩邊進行積分得到

(13)

將各變量用方波脈沖級數(shù)展開并代入式(13)，根據(jù)方波脈沖函數(shù)的脫關(guān)性，可得

(14)

由于

(15)

因而有

(16)

由于上述方程對t∈[0,T)內(nèi)的任何t值均成立，所以令等號兩邊對應(yīng)系數(shù)相等，可得

(17)

(18)

故有

(19)

式中

對每個不同的初始狀態(tài)向量可以求出N個子區(qū)間內(nèi)不同的狀態(tài)響應(yīng)，由于X0為n×n滿秩矩陣，可得X·i也是n×n滿秩矩陣，由式(19)可以解出A1為

(20)

由于

(21)

因而可得系統(tǒng)矩陣Ai的遞推公式

(22)

1.3 周期時變振動系統(tǒng)剛度與阻尼矩陣識別

p自由度阻尼、剛度周期時變的振動系統(tǒng)動力微分方程為

(23)

系統(tǒng)在做自由振動時，將式(23)用狀態(tài)方程來描述為

(24)

(25)

(26)

由式(25)可見，A中后p行為待求參數(shù)，為減少計算量，提高求解速度，將矩陣進行分塊

從而根據(jù)式(22)可推導(dǎo)出剛度、阻尼參數(shù)識別遞推公式(27)～(30)

(27)

(28)

(29)

(30)

對于具有周期時變參數(shù)的振動系統(tǒng)，根據(jù)Sylvester與Fourier理論，其自由響應(yīng)的頻譜為

(31)

2 二自由度周期時變仿真系統(tǒng)參數(shù)識別

圖1所示一個二自由度質(zhì)量-彈簧-阻尼時變系統(tǒng)，質(zhì)量m1=1 kg，m2=2 kg，剛度k1=1×104N/m，kc=2×104N/m，kv=5×103N/m，剛度時變頻率Ω=20 rad/s、初相位φ=π/2。時變剛度k2=kc+kvcos(Ωt+φ)。為了突出研究剛度變化，現(xiàn)假設(shè)阻尼c1=c2=0，外界激振力f1=f2=0。

圖1 二自由度剛度周期時變振動系統(tǒng)

圖2 二自由度剛度周期時變系統(tǒng)振動位移響應(yīng)

對x2位移信號進行頻譜分析，如圖3所示。

圖3 x2位移頻譜圖

本文共進行四組不同初始狀態(tài)(且線性無關(guān))下的響應(yīng)計算，根據(jù)遞推公式(27)～(30)，利用MATLAB語言編制參數(shù)識別程序，時變剛度識別結(jié)果及相對誤差如圖4～圖6所示。

圖4 時變剛度k2識別結(jié)果(時域圖)

圖5 時變剛度k2識別結(jié)果(頻域圖)

圖6 剛度識別相對誤差

為研究不同信噪比、計算步長對識別結(jié)果的精度影響，定義識別結(jié)果的平均絕對百分誤差(MAPE)為

(32)

信噪比與計算步長對識別精度的影響分別如表1、表2所示。

通過表1、表2可知，信號噪聲對識別結(jié)果有一定的影響，故在信號采集過程中應(yīng)盡量提高信噪比；計算步長越小，識別效果越好，但相應(yīng)的計算量也增大，故應(yīng)根據(jù)實際要求選擇合適的步長。

3 非對稱轉(zhuǎn)子剛度識別

為了對本文提出的參數(shù)識別方法進行實驗驗證，搭建具有剛度周期時變特性的非對稱(方軸)轉(zhuǎn)子實驗系統(tǒng)，如圖7所示；實驗所用儀器設(shè)備如表3所示。對非對稱轉(zhuǎn)子系統(tǒng)進行錘擊實驗(用力錘與水平成45°敲擊軸中心，即盤旁)，利用LMS Test.lab軟件進行數(shù)據(jù)采集與處理。

圖7 非對稱轉(zhuǎn)子實驗系統(tǒng)

非對稱轉(zhuǎn)子結(jié)構(gòu)如圖8、圖9所示，轉(zhuǎn)軸材料為30CrMnSi，輪盤材料為45#鋼，相關(guān)參數(shù)為：l=475 mm，d軸=14 mm，h=8 mm，ρ軸=7 750 kg/m3，E軸=2.01×1011Pa，μ軸=0.227，D盤=212 mm，ρ盤=7 850 kg/m3。為了實現(xiàn)軸的非對稱性，在軸上進行部分切割處理，使其在相互垂直的兩個方向上抗彎剛度不同(方軸)。實驗中，調(diào)節(jié)轉(zhuǎn)速控制器，使轉(zhuǎn)子在n=720 r/min下穩(wěn)定運轉(zhuǎn)，這樣軸在水平、垂直方向上的抗彎剛度則隨轉(zhuǎn)速時變；用力錘對轉(zhuǎn)子系統(tǒng)進行沖擊響應(yīng)實驗，通過降低轉(zhuǎn)子不平衡度與濾波的方法來消除轉(zhuǎn)子在離心力下的不平衡響應(yīng)影響，突出轉(zhuǎn)子在脈沖激勵后的自由響應(yīng)。

圖8 非對稱轉(zhuǎn)子結(jié)構(gòu)示意圖

圖9 非對稱轉(zhuǎn)子實物與測點分布圖

對非對稱轉(zhuǎn)子系統(tǒng)進行多次錘擊測試，利用電渦流位移傳感器采集軸中心的振動位移(水平、豎直方向)，用激光速度測試儀采集軸中心水平方向振動速度，用加速度傳感器采集轉(zhuǎn)子支承處的振動加速度(水平、豎直、軸向)。某次測試水平方向位移與速度響應(yīng)分別如圖10、圖11所示。

圖10 測點8水平方向位移時域響應(yīng)

圖11 測點9水平方向速度時域響應(yīng)

利用本文的方法對實驗中采集到的信號進行分析，從而得到非對稱轉(zhuǎn)子的時變剛度，如圖12所示。

圖12 非對稱轉(zhuǎn)子水平方向剛度變化曲線

可見運用本文方法識別出的非對稱轉(zhuǎn)子水平剛度由常值和周期時變兩部分組成，由于數(shù)據(jù)采集中噪聲的存在，識別出的剛度曲線存在一些毛刺。對該曲線進行頻域分析，如圖13所示。

由圖13可見，本文所識別出的剛度為與kc=1.25×105N/m，kv=3.76×104N/m，與利用ANSYS有限元軟件計算出的kc=1.27×105N/m，kv=3.82×104N/m相近。非對稱轉(zhuǎn)子水平剛度變化具有明顯的周期性，頻率成分為f=24 Hz，恰為轉(zhuǎn)子轉(zhuǎn)頻fn=720/60=12的二倍，與轉(zhuǎn)子轉(zhuǎn)一周剛度時變二次相符。不同計算步長(即不同采樣頻率)對剛度識別MAPE的影響如表4所示。可見，當選擇合適的步長(本實驗選取10-4s)時，平均絕對百分誤差小于1%。

圖13 非對稱轉(zhuǎn)子水平方向剛度頻譜圖

4 結(jié) 論

(1)本文所推導(dǎo)出的周期時變轉(zhuǎn)子系統(tǒng)參數(shù)識別遞推公式具有結(jié)構(gòu)簡明、計算省時等優(yōu)點。

(2)信號噪聲對識別結(jié)果有一定的影響，故在信號采集過程中應(yīng)盡量提高信噪比；計算步長越小，識別效果越好，但相應(yīng)的計算量也增大，故應(yīng)根據(jù)實際要求選擇合適的步長。

(3)通過對二自由度數(shù)值仿真模型和實際周期時變轉(zhuǎn)子系統(tǒng)分別進行參數(shù)識別，均較準確的識別出了時變剛度參數(shù)，其平均絕對百分誤差分別小于0.5%與1%，從而驗證了本方法的正確性和有效性。

(4)本文為齒輪、裂紋轉(zhuǎn)子等具有周期時變參數(shù)的工程機械結(jié)構(gòu)的參數(shù)識別及故障診斷提供一種新方法，有較大的工程應(yīng)用價值。