含間隙非線性二元翼段的系統(tǒng)辨識

孫玉凱，楊超，吳志剛

（北京航空航天大學 航空科學與工程學院，北京100083）

間隙非線性是一類重要的結構非線性，主要存在于具有相對運動的不同結構部件連接處［1］。實際結構在加工裝配過程中可能由于誤差和公差，其接觸表面往往會存在間隙；另外，在運行中結構的失效或松動也可能產(chǎn)生間隙。間隙非線性的存在使得結構表現(xiàn)出豐富的非線性動力學特性，如超臨界Hopf分岔和亞臨界Hopf分岔［2］、極限環(huán)振蕩及混沌現(xiàn)象［3］。在氣動彈性系統(tǒng)中，控制面間隙非線性的出現(xiàn)會在低于線性顫振邊界內(nèi)引起結構多種模式的極限環(huán)振蕩［4］，這種極限環(huán)振蕩現(xiàn)象將加速結構失效。由于實際結構中間隙大小不易測量，通常需要通過系統(tǒng)辨識來建立結構的數(shù)學模型［5］。

近年來，非線性動力學系統(tǒng)的系統(tǒng)辨識的發(fā)展迅速從學術研究向實際工程應用推廣，新發(fā)展的多自由度非線性系統(tǒng)辨識算法面向復雜的非線性結構。在動力學系統(tǒng)的非線性辨識領域，Kerschen等［6］總結了近年來的成果與發(fā)展，指出該領域未來的研究重點。對于間隙非線性動力學系統(tǒng)的非線性系統(tǒng)辨識，雖然國內(nèi)外學者做了大量的研究［7-10］，但缺乏準確估計間隙邊界的方法，不同的間隙辨識模型的結果差異較大；另外，間隙的存在將非線性動力學系統(tǒng)分成了多個分段線性系統(tǒng)，現(xiàn)有方法難以準確獲得每個分段線性系統(tǒng)的特性。

在工程應用中，通過地面振動試驗獲取的模態(tài)參數(shù)是結構描述和動力學分析的基礎。然而，在實際結構中，某一位置間隙的存在使結構整體表現(xiàn)出非線性特性，使地面振動試驗的測試結果產(chǎn)生偏差，無法獲得結構標稱線性系統(tǒng)特性。雖然諸如力-狀態(tài)圖法［11］等方法可以直觀地描述間隙恢復力狀態(tài)，但獲取標稱線性系統(tǒng)的模態(tài)特性十分困難。因此，能夠同時辨識標稱線性系統(tǒng)和非線性系數(shù)的非線性系統(tǒng)辨識方法則擁有巨大的工程應用前景。

本文針對含間隙非線性的動力學系統(tǒng)，利用條件逆譜（Conditioned Reverse Path，CRP）法和時域非線性子空間（Time-domain Nonlinear Subspace Identification，TNSI）法，分別建立了間隙非線性的數(shù)學描述，并成功辨識出非線性結構的標稱線性系統(tǒng)與非線性系數(shù)。通過對含間隙非線性的二元翼段進行系統(tǒng)辨識，驗證并對比了2種辨識方法。

1 非線性系統(tǒng)辨識方法

1.1 條件逆譜法

條件逆譜法［12］是對逆譜法［13］的推廣，是一種頻域系統(tǒng)辨識方法。在頻域內(nèi)從測量信號中分離出標稱線性系統(tǒng)的成分，從而獲得非線性系統(tǒng)的標稱線性系統(tǒng)的頻響函數(shù)（Frequency Response Function，F(xiàn)RF）與非線性系數(shù)。

對于一個包含N個自由度的非線性系統(tǒng)，其動力學方程可以表示為［12］

式中：M、C和K分別為系統(tǒng)的質量、阻尼和剛度矩陣；x（t）為位移/廣義位移向量；f（t）為力/廣義力向量；假設該系統(tǒng)包含n個非線性項，其中第j個非線性項的數(shù)學模型為yj（t），Lj為非線性位置向量，表示第j個非線性項所在的自由度，μj為該非線性項系數(shù)。

將式（1）通過傅里葉變換轉換到頻域：

式中：X（ω）、Yj（ω）、F（ω）分別為x（t）、yj（t）、f（t）的傅里葉變換；B（ω）=-ω2M+iωC+K為該系統(tǒng)標稱線性系統(tǒng)的動剛度矩陣。

傳統(tǒng)的模態(tài)參數(shù)辨識方法利用系統(tǒng)的頻響函數(shù)矩陣H（ω）=B（ω）-1辨識得到系統(tǒng)的模態(tài)參數(shù)。

對于多自由度系統(tǒng)，頻響函數(shù)矩陣通常利用單邊功率譜密度估計，即H1估計和H2估計，分別為

式中：G為單邊功率譜密度矩陣，下標分別表示對應的輸入、輸出通道，其中F表示力信號，X表示位移信號，如GFX表示力與位移的互功率譜。

條件逆譜法將F（ω）作為系統(tǒng)的輸出，將X（ω）和Yj（ω）作為系統(tǒng)的輸入，則式（2）可以轉換為

系統(tǒng)如圖1所示。

考慮到通過測量得到的x（t）和f（t）及通過響應計算得到的非線性響應yj（t），其對應的傅里葉變換信號中各非線性項之間互相關聯(lián)。以第2項非線性項的傅里葉變換Y2為例，其中包含與第1項 非 線 性 項 Y1相 關 的 Y2（+1）和 不 相 關 的Y2（-1）。由于存在 相關性，則Y1與Y2（+1）之間的關系可以由傳遞函數(shù)T12來表示，即

關系如圖2所示。

圖1 條件逆譜法系統(tǒng)示意圖Fig.1 Schematic diagram of system for CRP method

圖2 非線性項Y2 的分解Fig.2 Decomposition of nonlinear function Y2

類似的，對于第3項非線性項Y3的分解為

式中：Y3（+1）和Y3（+2）分別為第3項非線性項與第1、2項非線性項的相關成分，其關系分別通過T13和T23來表示。

一般情況下，第i項非線性項Yi中與前i-1個非線性項均無關的成分可以展開為

式中：類似式（3），Tji可以估計為

于是，位移響應信號X和激勵信號F均可通過上述方法展開，如圖3和圖4所示。

經(jīng)過分解得到的位移響應X（-1∶n）與所有非線性項無關，即可看作是標稱線性系統(tǒng)的位移響應為

圖3 位移響應信號分解Fig.3 Decomposition of displacement response signal

圖4 激勵信號分解Fig.4 Decomposition of excitation signal

同理，分解后的激勵響應F（-1∶n）為

1.1.1 標稱線性系統(tǒng)頻響函數(shù)估計

利 用 激 勵 信 號 F（-1∶n）和 位 移 響 應 信 號X（-1∶n）的單邊功率譜密度，參考H1和H2估計方法，可以估計標稱線性系統(tǒng)的頻響函數(shù)矩陣。

在進行標稱線性系統(tǒng)頻響函數(shù)估計的過程中，需要先計算非條件功率譜密度，并且需要注意的是，根據(jù)互惠條件，Gij=GHji，上標H表示矩陣的Herm itian變換。一般情況下，系統(tǒng)的條件功率譜密度為（r＜i，j）

式中：T為采樣周期。

當計算條件功率譜密度矩陣時，下標i和j可以替換成X或F計算?；跅l件逆譜法估計的頻響函數(shù)，可以得到該標稱線性系統(tǒng)的模態(tài)參數(shù)。

1.1.2 非線性系數(shù)辨識

根據(jù)式（14）可知，非線性系數(shù)使用遞歸方法從（μnLn）開始計算至（μ1L1）。

1.1.3 累積相干函數(shù)

在使用條件逆譜法的過程中，由于需要將非線性項看作系統(tǒng)的輸入，需要事先給定非線性項yj（t）的形式。為評估選擇的yj（t）描述是否符合系統(tǒng)實際，為此引入累積相干函數(shù)［14］作為評價標準。累積相干函數(shù)i為

1.2 時域非線性子空間法

時域非線性子空間法由Marchesiello和Garibaldi［15-16］總結并引入到非線性系統(tǒng)辨識中，該方法的關鍵在于將外力作為待辨識系統(tǒng)的輸入，其原理如圖5所示。子空間算法是通過諸如正交三角（QR）分解和奇異值分解（SVD）等矩陣運算辨識得到系統(tǒng)的狀態(tài)空間模型，屬于魯棒性強的非迭代算法。

圖5 非線性振動系統(tǒng)的反饋表示Fig.5 Feedback interpretation of nonlinear vibration systems

將式（1）中的非線性項移動到等式右邊，即

將式（16）改寫為一階狀態(tài)空間矩陣的形式，假設系統(tǒng)的狀態(tài)變量為z=［xT˙xT］T，則系統(tǒng)動力學方程的狀態(tài)空間表示為

式中：u（t）=［fTy1… yn］T表示系統(tǒng)的擴展輸入矩陣；下標c指代連續(xù)狀態(tài)空間系統(tǒng)，Ac、Bc、Cc、Dc分別為系統(tǒng)的狀態(tài)矩陣、輸入矩陣、輸出矩陣、直接傳遞矩陣，即

對上述連續(xù)系統(tǒng)離散采樣，假設采樣率為fs=1/T，則離散后的狀態(tài)空間方程為

式中：下標d表示離散系統(tǒng)。

根據(jù)數(shù)字采樣的零階保持假設，連續(xù)域與離散域之間的關系為

根據(jù)式（18），定義非線性系統(tǒng)擴展的頻響函數(shù)矩陣HE為［16］

考慮

將式（22）代入到式（21）中，有

式中：P12=（K+jωC-ω2M）-1M。

系統(tǒng)的標稱線性系統(tǒng)的頻響函數(shù)為

于是，式（23）最終簡化為

根據(jù)式（25）可以辨識得到標稱線性系統(tǒng)的頻響函數(shù)及非線性系數(shù)。

2 數(shù)值算例

在空氣動力學與氣動彈性等相關領域中，二元翼段通常是驗證新方法、檢驗新技術手段的典型模型，其結構自由度較低，并且可使用二維氣動力開展研究。對于結構間隙引起的非線性氣動彈性問題，含間隙二元翼段也是十分理想的研究對象，不僅有解析解或半解析解作為對照［17］，更得到了國內(nèi)外學者的廣泛研究［18］；在風洞試驗的研究中，含間隙的非線性翼段模型也受到了廣泛的關注［19］。獲取準確的間隙參數(shù)是理論分析與試驗探究相結合的關鍵環(huán)節(jié)，因此，本文將以含間隙的二元翼段為例，開展非線性系統(tǒng)辨識的研究，討論對比2種非線性辨識方法。

2.1 非線性二元翼段模型

考慮一個具有沉浮h和俯仰α兩個自由度的二元翼段［20］，如圖6所示。圖中：b為二元翼段半弦長，ab為剛心E到翼段中心O的距離（剛心在中心后為正），xa為重心G到剛心E的距離（重心在剛心后為正）。h為沉浮自由度位移，α為俯仰自由度位移。在俯仰自由度上存在間隙非線性，忽略系統(tǒng)的阻尼。沉浮自由度和俯仰自由度的剛度分別為kh和kα。

圖6 二自由度二元翼段模型Fig.6 2-DOF wing section model

通過拉格朗日方程可以寫出其動力學方程為式中：m為翼段質量；Sα為翼段對剛心E軸的靜矩；Iα為翼段對剛心E軸的轉動慣量；kh為沉浮自由度的剛度系數(shù)；fh和fα分別表示對應自由度上的外力；fnl（α）為包含間隙的俯仰自由度非線性恢復力，其數(shù)學描述為

其中：kα為俯仰自由度剛度系數(shù)；dα為間隙大小。式（27）表示的間隙模型如圖7所示。

圖7 間隙非線性的恢復力曲線Fig.7 Restoring force of freeplay nonlinearity

2.2 標稱線性系統(tǒng)

由于間隙的存在，式（26）描述的非線性系統(tǒng)可以拆分成不同的線性系統(tǒng)的組合。不妨定義該間隙非線性系統(tǒng)的2個標稱線性系統(tǒng)分別是上標稱線性系統(tǒng)（Overlying Linear System，OLS）和下標 稱 線 性 系 統(tǒng) （Underlying Linear System，ULS）［21］，其定義分別如下：

1）OLS。當非線性間隙不存在時，即式（27）中的dα=0，此時的系統(tǒng)稱為上標稱線性系統(tǒng)，即普遍認為的間隙非線性系統(tǒng)中的標稱線性系統(tǒng)。式（26）的OLS為

2）ULS。當kα=0時，即式（27）描述的非線性恢復力fnl（α）=0，此時的系統(tǒng)稱為下標稱線性系統(tǒng)。式（26）的ULS為

該二元翼段的結構基本參數(shù)如表1所示，其沉浮自由度剛度為kh=782.74 N/m，俯仰自由度剛度為kα=24.22 N·m/rad，若俯仰自由度間隙dα=0 rad，即得到該二元翼段的標稱線性情況。該二元翼段標稱線性情況的固有頻率為沉浮固有頻率2.61 Hz，俯仰固有頻率為5.05 Hz。

2.3 數(shù)值仿真與模態(tài)辨識

對該二元翼段的俯仰自由度施加頻帶范圍0～50 Hz的隨機激勵，隨機激勵的強度由激勵信號的均方根（Root Mean Square，RMS）表示，使用Runge-Kutta數(shù)值仿真得到二元翼段的振動響應。仿真采樣率設置為1 000 Hz，仿真總時長為120 s。非線性系統(tǒng)仿真時的間隙大小取dα=0.01 rad。為模擬真實的試驗情況，仿真中添加2%的噪聲。取響應穩(wěn)定之后的一段長度為214個樣本點的數(shù)據(jù)，利用平均周期圖法及根據(jù)式（3）估計該二元翼段的頻響函數(shù)曲線。

在不同的激勵水平作用下（激勵水平對應的隨機力RMS值見表2），系統(tǒng)的頻響函數(shù)曲線如圖8所示。在選取的0～15 Hz的頻率范圍內(nèi)，隨著激勵的增加，頻響函數(shù)曲線呈現(xiàn)出2個穩(wěn)定的峰值，分別對應二元翼段的沉浮模態(tài)和俯仰模態(tài)。其中，由于俯仰自由度間隙非線性的影響，使得頻響函數(shù)曲線中俯仰模態(tài)的峰值隨激勵水平變化（見圖8中藍色箭頭）。

表2 激勵水平與對應的隨機力RM S值Tab le 2 Excitation levels and corresponding RM S value of random excitation force

圖9顯示了頻帶0～30 Hz內(nèi)，低激勵水平（隨機力RMS為0.0048N·m）與高激勵水平（隨機力RMS為0.230 8 N·m）的頻響函數(shù)曲線對比。低激勵幅值的頻響函數(shù)曲線畸變嚴重，無法準確獲得系統(tǒng)的模態(tài)特性；而且由于間隙非線性的存在，使得系統(tǒng)響應出現(xiàn)高頻振動。

圖8 隨激勵水平變化的頻響函數(shù)Fig.8 Variation of FRF with excitation level

3 非線性系統(tǒng)辨識結果

3.1 條件逆譜法辨識結果

在使用條件逆譜法辨識間隙非線性時，非線性函數(shù)fnl（α）的描述定義為

式中：μ為待辨識的非線性系數(shù)；間隙非線性模型通過參數(shù)ν來調(diào)整，該數(shù)學模型描述的準確性通過累積相干函數(shù)評估。

為了通過系統(tǒng)辨識獲得非線性系統(tǒng)的OLS，根據(jù)式（30）將描述非線性的函數(shù)y定義為

當激勵水平為隨機力RMS 0.048 6 N·m時，利用條件逆譜法辨識獲得系統(tǒng)的頻響函數(shù)曲線如圖10所示?？梢钥闯?，條件逆譜法可以獲得清晰的OLS頻響函數(shù)曲線。根據(jù)條件逆譜法獲得的頻響函數(shù)曲線，利用模態(tài)參數(shù)辨識方法即可獲得標稱線性系統(tǒng)的模態(tài)參數(shù)。

利用牛頓切線法對非線性描述式（31）中的ν進行迭代優(yōu)化，使得辨識的累積相干函數(shù)最大。俯仰自由度下的累積相干函數(shù)如圖11所示，其均值為0.95，表明辨識得到的非線性模型與實際非線性系統(tǒng)特性吻合，此時ν=261。

在得到標稱線性系統(tǒng)的頻響函數(shù)矩陣之后，利用式（14）可以得到系統(tǒng)的非線性系數(shù)，辨識得到的非線性系數(shù)μ的結果如圖12所示。表3給出了在不同激勵水平下辨識得到的非線性系數(shù)。可以看出，非線性系數(shù)的辨識結果一致性較好。通過式（31）重構非線性恢復力曲線如圖13所示?？梢钥闯?，辨識得到的非線性恢復力曲線與真實的恢復力曲線一致。

圖10 條件逆譜法辨識得到的OLS頻響函數(shù)曲線Fig.10 FRF of OLS identified by CRPmethod

圖11 俯仰自由度的累積相干函數(shù)Fig.11 Cumulative coherence function of pitch DOF

圖12 條件逆譜法辨識的非線性系數(shù)Fig.12 Nonlinear coefficient identified by CRPmethod

表3 不同激勵水平下辨識的非線性系數(shù)Tab le 3 Identified non linear coefficient under different excitation levels

圖13 條件逆譜法辨識的非線性恢復力曲線Fig.13 Nonlinear restoring force identified by CRPmethod

3.2 時域非線性子空間法辨識結果

在使用時域非線性子空間法辨識間隙非線性時，非線性函數(shù)fnl（α）的定義參考式（27），選定20個非線性項，定義如下：

式中：dj為可選定的間隙大??；μj為第j個非線性項的非線性系數(shù)。

為了能夠真實地刻畫間隙大小，將測量的響應位移范圍［0，max（|α|）］均分為20段，分別為［dj-1，dj］，d0=0，j=1，2，…，20［16］。

當激勵水平為隨機力RMS 0.022 6N·m時，利用時域非線性子空間法辨識得到的非線性系數(shù)μj如圖14所示?？梢钥闯?，非線性系數(shù)主要集中在0.01 rad附近。

為了得到更加精確的結果，可以進一步縮小選定的間隙dj的區(qū)間為［0.008，0.012］rad，將該區(qū)間均分為10段，通過第二次時域非線性子空間法辨識得到的非線性系數(shù)如圖15所示。圖16給出了上述2次辨識獲得的ULS的頻響函數(shù)曲線，圖17顯示了2次辨識獲得的非線性恢復力曲線與標稱曲線的對比。結果表明，時域非線性子空間辨識可以準確地獲得間隙參數(shù)，對非線性恢復力曲線中拐點的估計也十分準確。

利用辨識得到的非線性二元翼段ULS與間隙參數(shù)，可以重構該非線性二元翼段的標稱線性系統(tǒng)，其頻響函數(shù)曲線如圖18所示。可以看出，時域非線性子空間法可以準確地辨識出含間隙非線性系統(tǒng)的標稱線性系統(tǒng)。子空間辨識法獲得的標稱線性系統(tǒng)的頻響函數(shù)可以方便地使用Poly-MAX［22］等模態(tài)參數(shù)辨識方法辨識，圖19的穩(wěn)態(tài)圖可以直觀地反映系統(tǒng)的模態(tài)特性。

圖14 區(qū)間［0，max（|α|）］內(nèi)的非線性系數(shù)分布Fig.14 Nonlinear parameter distribution in theinterval［0，max（|α|）］

圖15 區(qū)間［0.008，0.012］rad內(nèi)的非線性系數(shù)分布Fig.15 Nonlinear parameter distribution in the interval［0.008，0.012］rad

圖16 辨識得到的非線性二元翼段ULS的頻響函數(shù)曲線Fig.16 Identified FRF of ULS of nonlinear 2-DOF wing section

圖17 時域非線性子空間法辨識的非線性恢復力曲線Fig.17 Nonlinear restoring force identified by TNSImethod

圖18 辨識得到的非線性二元翼段OLS的頻響函數(shù)曲線Fig.18 Identified FRF of OLS of nonlinear 2-DOF wing section

圖19 辨識得到的OLS頻響函數(shù)曲線穩(wěn)態(tài)圖Fig.19 Stabilization diagram of identified FRF of OLS

3.3 方法對比

條件逆譜法和時域非線性子空間法都能準確地辨識出非線性系統(tǒng)的OLS。一方面，由于2種方法選用的非線性數(shù)學描述不同，條件逆譜法可以直接獲得系統(tǒng)的OLS，但對間隙非線性拐點的辨識不太準確；時域非線性子空間法可直接獲得系統(tǒng)的ULS，也能準確地辨識間隙非線性拐點，系統(tǒng)的OLS則需要重構。另一方面，條件逆譜法辨識獲得的非線性系數(shù)隨頻率線有波動，而由于子空間辨識算法的魯棒性較好，在非線性特性較強時也能辨識出清晰的頻響函數(shù)曲線。2種方法的特性在表4中做了系統(tǒng)的對比和總結。

表4 條件逆譜法和時域非線性子空間法屬性對比Tab le 4 Proper ty com parison between CRP and TNSI

4 結 論

1）條件逆譜法和時域非線性子空間法均可辨識含間隙非線性動力學系統(tǒng)，獲得系統(tǒng)的標稱線性系統(tǒng)與對應的非線性系數(shù)。

2）2種方法辨識過程中對間隙非線性采用了不同的數(shù)學描述，均能獲得與理論值一致的非線性恢復力曲線。時域非線性子空間法對間隙非線性拐點的辨識更加準確，條件逆譜法可以直接獲得標稱線性系統(tǒng)特性。