張 波 張文博 彭志科 ,?,

* (寧夏大學(xué)機械工程學(xué)院,銀川 750021)

? (上海交通大學(xué)機械系統(tǒng)與振動國家重點實驗室,上海 200240)

引言

頻率響應(yīng)函數(shù)是刻畫線性系統(tǒng)的一個基本概念,它是經(jīng)典控制論和系統(tǒng)設(shè)計的理論基礎(chǔ),在工程實踐中已得到了廣泛應(yīng)用[1].考慮到利用頻率響應(yīng)函數(shù)對線性系統(tǒng)分析帶來的獨特優(yōu)勢,一些學(xué)者基于Volterra 級數(shù)理論,將線性系統(tǒng)中頻率響應(yīng)函數(shù)的概念拓展到了非線性系統(tǒng)中,其中廣義頻率響應(yīng)函數(shù)(generalized frequency response functions,GFRFs)是非線性系統(tǒng)中各階Volterra 核函數(shù)的多維傅里葉變換[2].廣義頻率響應(yīng)函數(shù)是描述非線性系統(tǒng)頻域特性的一系列多維函數(shù),非線性系統(tǒng)的第一階廣義響應(yīng)函數(shù)與其線性派生系統(tǒng)的頻率響應(yīng)函數(shù)相同.從20 世紀(jì)70 年代初以來,廣義頻率響應(yīng)函數(shù)已經(jīng)應(yīng)用于非線性系統(tǒng)的頻率特性研究,如頻率響應(yīng)的諧波特性、增益壓縮擴張?zhí)匦院拖嗷フ{(diào)制特性等[3];在機械工程和生物工程等領(lǐng)域中廣義頻率響應(yīng)函數(shù)也用于分析系統(tǒng)的非線性特性[4-6].由于廣義頻率響應(yīng)函數(shù)的多維特性,導(dǎo)致高階廣義頻率響應(yīng)函數(shù)物理意義難以理解,同時高階廣義頻率響應(yīng)函數(shù)也不易被直觀表示出來,這些因素限制了廣義頻率響應(yīng)函數(shù)在實際工程應(yīng)用中的發(fā)展.為了解決這個問題,在過去的幾十年里,人們提出了幾種非線性系統(tǒng)的一維頻域表示方法.例如Lang 等[7-8]提出的非線性輸出頻率響應(yīng)函數(shù)(nonlinear output frequency response functions,NOFRFs)、非線性系統(tǒng)的輸出頻率響應(yīng)函數(shù)(output frequency response function,OFRF)以及Rijlaarsdam 等[9]提出的高階正弦輸入表示函數(shù)(higher order sinusoidal input describing functions,HOSIDF).已有研究表明,輸出頻率響應(yīng)函數(shù)用于在頻域中對非線性系統(tǒng)的設(shè)計[10-11],高階正弦輸入描述函數(shù)用于研究受正弦輸入影響的非線性系統(tǒng)[12].而非線性輸出頻率響應(yīng)函數(shù)則是一系列關(guān)于頻率的一維函數(shù),它對非線性系統(tǒng)頻域分析提供了一系列類似于伯德(Bode)圖的頻域表示方法.此外,基于非線性輸出頻率響應(yīng)函數(shù)的非線性系統(tǒng)分析包括對所研究的非線性系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)進行計算,用非線性輸出頻率響應(yīng)函數(shù)作為相關(guān)指標(biāo)來揭示系統(tǒng)的動態(tài)特性[13-14].目前,非線性輸出頻率響應(yīng)函數(shù)已經(jīng)成功應(yīng)用于工程結(jié)構(gòu)系統(tǒng)的狀態(tài)監(jiān)測和故障診斷[15-24].

研究表明,利用非線性輸出頻率響應(yīng)函數(shù)進行非線性系統(tǒng)分析及結(jié)構(gòu)損傷檢測的關(guān)鍵是如何準(zhǔn)確計算出系統(tǒng)的各階非線性輸出頻率響應(yīng)函數(shù).目前,有關(guān)非線性輸出頻率響應(yīng)函數(shù)的計算方法大多數(shù)是基于系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)對輸入激勵幅值變化不敏感性的原理進行計算,即對系統(tǒng)進行多次相同頻率不同幅值的激勵,得到系統(tǒng)輸出響應(yīng)后利用最小二乘法(least squares method,LSM)計算出系統(tǒng)的各階非線性輸出頻率響應(yīng)函數(shù).而在用最小二乘方法計算非線性系統(tǒng)的高階非線性輸出頻率響應(yīng)函數(shù)時,往往需要輸入激勵次數(shù)不小于系統(tǒng)的截斷階數(shù),而且在結(jié)構(gòu)損傷檢測中大多采用離線檢測的方式.此外,Feijoo 等[25]提出了非線性系統(tǒng)伴隨線性方程(associated linear equations,ALE)的概念,Bayma 等[26]根據(jù)伴隨線性方程提出了基于帶外部輸入的非線性自回歸模型(nonlinear auto regressive with eXogenous input,NARX)非線性輸出頻率響應(yīng)函數(shù)的計算方法,但該方法計算過程中涉及丟番圖方程的求解問題,同時由于伴隨線性方程法僅適用于非線性項為輸入激勵或輸出響應(yīng)多項式的非線性系統(tǒng),所以無法求解含有輸入激勵或輸出響應(yīng)導(dǎo)數(shù)項的非線性系統(tǒng),如范德波爾系統(tǒng)等.最近,Zhu等[27]提出了基于NARX 模型表示的廣義伴隨線性方程法(generalized associated linear equations,GALEs),并且將該方法應(yīng)用于機床刀具損傷狀態(tài)監(jiān)測及火車車輪疲勞損傷的在線檢測.

從非線性系統(tǒng)頻域分析的研究情況來看:關(guān)于非線性輸出頻率響應(yīng)函數(shù)的計算方法大多數(shù)基于最小二乘法.雖然伴隨線性方程法可以準(zhǔn)確地計算出非線性系統(tǒng)任意階非線性輸出響應(yīng),但其適用范圍卻有很大局限性,所以進一步研究伴隨線性方程是利用非線性輸出頻率響應(yīng)函數(shù)進行非線性系統(tǒng)分析與結(jié)構(gòu)損傷檢測的基礎(chǔ).本文基于非線性微分方程(nonlinear differential equation,NDE)模型表示的非線性系統(tǒng)廣義頻率響應(yīng)函數(shù)遞歸計算公式及Volterra級數(shù)理論,對由NDE 模型表示的一類非線性系統(tǒng)的廣義伴隨線性方程進行研究,采用數(shù)值計算的方法對這類非線性系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)進行分析.同時,利用線性算子理論對非線性系統(tǒng)中典型非線性效應(yīng)的產(chǎn)生機理進行研究,為利用非線性輸出頻率響應(yīng)函數(shù)進行結(jié)構(gòu)損傷檢測及非線性系統(tǒng)的分析和設(shè)計提供了一種有效途徑.

1 非線性系統(tǒng)NDE 模型的廣義伴隨線性方程

1.1 廣義頻率響應(yīng)函數(shù)

根據(jù)Weierstrass 函數(shù)逼近理論[28],在封閉且有界區(qū)間上,任何連續(xù)函數(shù)都能用一組多項式函數(shù)來對其進行任意精度的一致逼近,所以工程實際中的許多非線性系統(tǒng)都可以描述為如下多項式型非線性系統(tǒng)[29-30]

式中,u(t) 為系統(tǒng)的輸入,y(t) 為系統(tǒng)的輸出,p+q=n,Cp,q(l1,l2,···,lp+q) 為多項式系數(shù),N為u(t) 與y(t)的最大非線性階數(shù),L為微分的最大階數(shù),求和符號為

微分算子 D 定義為

針對由式(1)所描述的系統(tǒng),當(dāng)系統(tǒng)輸入是具有實際物理意義的信號u(t) 時,該系統(tǒng)可產(chǎn)生確定的輸出響應(yīng)y(t),且系統(tǒng)必定存在非零線性輸出系數(shù),即C1,0(0)≠0 .若由NDE 模型表示單輸入單輸出非線性系統(tǒng),則式(1)可以表示為

式中,等式右邊Cp,q(l1,l2,···,lp+q) 滿足條件:C1,0(0)?Cp,q(l1,l2,···,lp+q) .

若式(4)表示的系統(tǒng)在零平衡點漸近穩(wěn)定,則該系統(tǒng)的輸出y(t) 可以用Volterra 級數(shù)表示為

式中,u(t) 為系統(tǒng)的輸入,y(t) 為系統(tǒng)的輸出,yn(t) 為系統(tǒng)的n階非線性輸出響應(yīng),hn(τ1,τ2,···,τn) 為非線性系統(tǒng)的n階廣義脈沖響應(yīng)函數(shù)或n階Volterra 核函數(shù),N為Volterra 級數(shù)的截斷階數(shù).對式(5)做傅里葉變換(Fourier transform,FT)后可得[31]

式中,Y(jω) 為系統(tǒng)輸出y(t) 的傅里葉變換,U(jω) 為系統(tǒng)輸入u(t) 的傅里葉變換,Yn(jω) 為系統(tǒng)的n階非線性輸出響應(yīng)yn(t) 的傅里葉變換.Hn(jω1,jω2,···,jωn)為系統(tǒng)n階Volterra 核函數(shù)hn(τ1,τ2,···,τn) 的n維傅里葉變換,被稱為廣義頻率響應(yīng)函數(shù).式(6) 中為被積函數(shù)在n維超平面 ω1+ω2+···+ωn=ω 上的多重積分.

Billings 等[29]根據(jù)諧波探測法推導(dǎo)出了由NDE 模型表示的計算非線性系統(tǒng)廣義頻率響應(yīng)函數(shù)遞歸表達(dá)式

式中,LN=(l1,l2,···,lp+q) 為與Cp,q(l1,l2,···,lp+q) 相對應(yīng)的系數(shù)指標(biāo).

1.2 非線性輸出頻率響應(yīng)函數(shù)

由式(6)可知非線性系統(tǒng)的廣義頻率響應(yīng)函數(shù)是多維函數(shù),各階廣義頻率響應(yīng)函數(shù)的維數(shù)等于它的階數(shù),由于廣義頻率響應(yīng)函數(shù)的階數(shù)超過3 階以后很難直觀表示出來,這對廣義頻率響應(yīng)函數(shù)在工程中的分析估計造成了極大困難.為了克服這個困難,Billings 等[7]提出了非線性系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)的概念.對于單輸入單輸出非線性系統(tǒng),非線性輸出頻率響應(yīng)函數(shù)為

式中,Yn(jω) 為系統(tǒng)的n階非線性輸出頻譜,Un(jω)為系統(tǒng)的n階非線性輸入頻譜,它是系統(tǒng)輸入u(t)的n次冪un(t) 的傅里葉變換,即

圖1 表示非線性系統(tǒng)中的各階非線性輸出頻率響應(yīng)函數(shù)與系統(tǒng)輸入輸出的關(guān)系,也揭示了線性系統(tǒng)中的頻率響應(yīng)函數(shù)與非線性系統(tǒng)中非線性輸出頻率響應(yīng)函數(shù)之間的關(guān)系.即n=N=1 時,Gn(jω)=G1(jω)表示線性系統(tǒng)的頻率響應(yīng)函數(shù).

圖1 非線性系統(tǒng)的NOFRFsFig.1 The NOFRFs of nonlinear systems

由式(9)描述的非線性輸出頻率響應(yīng)函數(shù)的最大特點是一維特性,它是廣義頻率響應(yīng)函數(shù)在不同加權(quán)因子作用下的加權(quán)平均值,所以非線性輸出頻率響應(yīng)函數(shù)對非線性系統(tǒng)在頻域的分析帶來了極大的便利,這也是非線性輸出頻率響應(yīng)函數(shù)在結(jié)構(gòu)狀態(tài)檢測及故障診斷等領(lǐng)域受到廣泛應(yīng)用的重要原因.

1.3 NDE 模型廣義伴隨線性方程

根據(jù)NDE 模型的廣義頻率響應(yīng)函數(shù)遞歸表達(dá)式(7)及非線性輸出頻率響應(yīng)函數(shù)的定義式(9),在式(7)左右兩邊同時乘并在n維超平面 ω1+ω2+···+ωn=ω 上積分得:

式(7)左邊結(jié)果為

式(7)右邊結(jié)果為

其中,An,u(jω) 為純輸入非線性項,且

式中An,uy(jω) 為輸入輸出耦合非線性項,且

式中

對式(17)兩邊同時做傅里葉逆變換得

An,y(jω) 為純輸出非線性項,且

式中,F T[·] 表示傅里葉變換.

分析可知,式(20)在表達(dá)形式上與線性系統(tǒng)的頻域表達(dá)形式類似,可以看作是該線性系統(tǒng)的頻率響應(yīng)函數(shù),Vn(jω) 是該線性系統(tǒng)廣義輸入的傅里葉變換,繼續(xù)對式(20)左右兩邊做傅里葉逆變換得

式中,Vn(t) 是Vn(jω) 的傅里葉逆變換,稱為廣義輸入,TL為線性算子.式(25)稱為非線性系統(tǒng)(4)的廣義伴隨線性方程[27].若在式(25)中已知線性算子 TL與Vn(t),則由NDE 模型描述的非線性系統(tǒng)(4) 的n階非線性輸出響應(yīng)yn(t) 可以通過求解簡單的線性方程得到.從廣義伴隨線性方程中計算出n階非線性輸出響應(yīng)yn(t) 后,可由式(11)計算出非線性系統(tǒng)的任意階非線性輸出頻率響應(yīng)函數(shù).

1.4 NDE 模型廣義伴隨線性方程的推導(dǎo)

對于由NDE 模型所描述的非線性系統(tǒng)(4),由廣義伴隨線性方程的表達(dá)式(13)～式(19)得

式中,n∈N?,LN=(l1,l2,···,lp+q) 是式(8)中與Cp,q(l1,l2,···,lp+q) 相對應(yīng)的指標(biāo).

式(26)稱為由NDE 模型所描述的非線性系統(tǒng)(4)的廣義伴隨線性方程,比較式(7)與式(26)可知兩者在表達(dá)形式上相似,并且式(27)中yn(t) 與可以理解為式(7) 中的廣義頻率響應(yīng)函數(shù)Hn(jω1,jω2,···,jωn) 與在時域中的類比,它揭示了一種關(guān)于非線性系統(tǒng)時域與頻域的新關(guān)系.

通過式(25)與式(26)可得,由NDE 模型表示的非線性系統(tǒng)廣義伴隨線性方程中的線性算子 TL與廣義輸出Vn(t) 表達(dá)式為

將式(27)代入式(29)可得Vn(t) 的表達(dá)式,其中廣義輸入Vn(t) 是由系統(tǒng)輸入激勵u(t) 及其各階微分Dliu(t) 和系統(tǒng)前n-1 階的輸出響應(yīng)yn-1(t),yn-2(t),···,y1(t) 及其各階微分 Dliyj(t) 的組合函數(shù).因此,通過求解式(26)表示的廣義伴隨線性方程后,可得系統(tǒng)(4)的各階非線性輸出響應(yīng)yn(t) .

1.5 廣義伴隨線性方程計算非線性輸出頻率響應(yīng)函數(shù)的步驟

廣義伴隨線性方程計算系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)分為以下3 個步驟:

(1) 將NDE 模型描述的非線性系統(tǒng)改寫成式(4)的形式;

(2) 由式(26)表示的廣義伴隨線性方程求解出系統(tǒng)的各階非線性輸出響應(yīng)yn(t) ;

(3) 由式(11)計算出系統(tǒng)的各階非線性輸出頻率響應(yīng)函數(shù)Gn(jω) .

1.6 非線性輸出頻率響應(yīng)函數(shù)的遞歸計算法與耦合計算法

計算NDE 模型的廣義伴隨線性方程(26)時,由式(25)可知,若廣義輸入Vn(t) 是由線性算子 TL的本征函數(shù)構(gòu)成的線性組合,則廣義伴隨線性方程(26)存在解析解.在實際工程系統(tǒng)中,系統(tǒng)的輸入激勵往往是沖擊、隨機或非平穩(wěn)信號,這類輸入激勵的廣義伴隨線性方程并不存在解析解,所以需要用數(shù)值計算的方法對廣義伴隨線性方程(26)進行求解.由1.5 節(jié)中廣義伴隨線性方程計算系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)的步驟可知,計算非線性輸出頻率響應(yīng)函數(shù)的關(guān)鍵是計算系統(tǒng)的廣義伴隨線性方程.

基于NDE 模型的廣義伴隨線性方程方法計算系統(tǒng)非線性輸出頻率響應(yīng)函數(shù)的一種方法是:在求解NDE 模型的廣義伴隨線性方程時,用數(shù)值計算(如Runge-Kutta)方法逐階求解廣義伴隨線性方程(26)可得系統(tǒng)的各階非線性輸出響應(yīng)yi(t),根據(jù)非線性輸出頻率響應(yīng)函數(shù)的定義式(11)可計算出系統(tǒng)前N階的非線性輸出頻率響應(yīng)函數(shù),該方法稱為遞歸計算法(recursive computational method,RCM),其計算步驟如表1 所示.

非線性輸出頻率響應(yīng)函數(shù)的另一種數(shù)值計算方法是:在求解NDE 模型的廣義伴隨線性方程時,將系統(tǒng)的前N階廣義伴隨線性方程視為具有坐標(biāo)耦合的N個線性方程組,用數(shù)值計算(如4 階Runge-Kutta)方法聯(lián)立求解耦合方程組后得到系統(tǒng)的各階非線性輸出響應(yīng)yi(t),最后根據(jù)非線性輸出頻率響應(yīng)函數(shù)的定義式(11)計算出系統(tǒng)前N階的非線性輸出頻率響應(yīng)函數(shù),該方法稱為耦合計算法(coupled computational method,CCM),其計算步驟如表2 所示.

表2 耦合計算法Table 2 Coupled computational method

2 數(shù)值仿真與分析

本節(jié)以范德波爾振子為例,利用數(shù)值仿真實驗對上述理論分析結(jié)果的正確性進行驗證.針對廣義伴隨線性方程求解非線性輸出頻率響應(yīng)函數(shù)的數(shù)值計算問題,分析了遞歸計算法與耦合計算法的計算結(jié)果產(chǎn)生差異的原因.最后,以達(dá)芬系統(tǒng)為例,通過廣義伴隨線性方程與線性算子理論分析了達(dá)芬系統(tǒng)中典型非線性效應(yīng)產(chǎn)生的原因.

2.1 范德波爾(van der Pol)系統(tǒng)

范德波爾振子的NDE 模型為[32]

式中,m,c,k為系統(tǒng)的線性參數(shù),CE為系統(tǒng)的非線性參數(shù),u(t) 為系統(tǒng)的輸入激勵,y(t) 為系統(tǒng)的輸出響應(yīng).

將式(30)代入式(4)得

將式(31)和式(32)代入式(26)得系統(tǒng)前7 階廣義伴隨線性方程為

設(shè)系統(tǒng)的參數(shù)為:m=1 kg,c=30 N·s/m,k=1.0×104N/m,CE=7.5×106N·s/m3,輸入激勵

式(30)的解y(t) 與式(33)中各階廣義伴隨線性方程的解yi(t) 均由4 階Runge-Kutta 方法計算得到,其中采樣頻率f s=1024 Hz .當(dāng)系統(tǒng)的最大非線性階數(shù)為N時,由式(5)可得系統(tǒng)輸出響應(yīng)yGN(t) 為

圖2 所示為Runge-Kutta 法與廣義伴隨線性方程法計算結(jié)果對比.圖2(a) 表示最大非線性階數(shù)N取不同值時系統(tǒng)頻域響應(yīng)對比.圖2(b)為N=7時,時域響應(yīng)對比.結(jié)果表明,取適當(dāng)?shù)慕財嚯A數(shù)N,廣義伴隨線性方程能夠很好地表示系統(tǒng)的輸出響應(yīng),從而驗證了由NDE 模型表示的非線性系統(tǒng)廣義伴隨線性方程(26)的有效性.

圖2 Runge-Kutta 與GALEs 響應(yīng)對比Fig.2 Comparison of the responses by using Runge-Kutta and GALEs

由式(33)可知,在零初始條件下,范德波爾振子式(30)的第2 階、4 階和6 階廣義伴隨線性方程的解y2(t),y4(t),y6(t) 為平凡解(0 解),即系統(tǒng)的前7 階非線性輸出頻率響應(yīng)函數(shù)只有G1(j2πf),G3(j2πf),G5(j2πf) 和G7(j2πf) 為非0 解.

圖3 表示用不同方法計算出的范德波爾系統(tǒng)前7 階非線性輸出頻率響應(yīng)函數(shù),其中GALEs-RCM 曲線和GALEs-CCM 曲線分別表示用遞歸計算法和耦合計算法的計算結(jié)果.為了比較不同方法計算非線性輸出頻率響應(yīng)函數(shù)的差異,這里與傳統(tǒng)的最小二乘方法計算非線性輸出頻率響應(yīng)函數(shù)作對比.由式(6)可知,非線性系統(tǒng)的輸出響應(yīng)頻譜為

圖3 范德波爾系統(tǒng)前7 階非線性輸出頻率響應(yīng)函數(shù)Fig.3 The first 7 nonlinear output frequency response functions of van der Pol system

最小二乘方法計算系統(tǒng)的前N階非線性輸出頻率響應(yīng)函數(shù)的表達(dá)式為[13]

其中,αm∈R+,表示系統(tǒng)輸入激勵的次數(shù),N表示系統(tǒng)的截斷階數(shù),.由式(33)可知范德波爾系統(tǒng)的前7 階非線性輸出響應(yīng)只有y1(t),y3(t),y5(t) 和y7(t) 為非零解,取n={1,3,5,7},αm={0.7,1,1.3,1.5},,N=4 .圖 3 中 的LSM 曲線為最小二乘法的計算結(jié)果.

圖3 表明最小二乘方法可以準(zhǔn)確地計算出低階非線性輸出頻率響應(yīng)函數(shù),而廣義伴隨線性方程法可以準(zhǔn)確地計算出任意高階非線性輸出頻率響應(yīng)函數(shù),但遞歸計算法對高階非線性輸出頻率響應(yīng)函數(shù)的計算結(jié)果偏差較大,這是因為用Runge-Kutta 方法逐階遞歸計算式(33)時,表1 中的第S2-1 步需要進行插值操作,目的是將離散化的yi(k) 進行連續(xù)化,便于求解高階的廣義伴隨線性方程,但同時也引入了插值偏差,低階廣義伴隨線性方程的數(shù)值計算偏差會隨著低階非線性輸出響應(yīng)一同作為高階廣義伴隨線性方程的輸入激勵,從而使高階廣義伴隨線性方程的數(shù)值計算偏差逐階累積,最終導(dǎo)致高階非線性輸出響應(yīng)計算結(jié)果偏差較大.由表2 中的第S2 步可知,耦合計算法在求廣義伴隨線性方程時并未進行函數(shù)插值操作,所以耦合計算法未引入插值誤差,比遞歸計算法的計算精度更高,最終非線性輸出頻率響應(yīng)函數(shù)的計算結(jié)果更準(zhǔn)確.

圖4 表示在耦合計算法中采用4 階Runge-Kutta 方法(ODE45)的求解結(jié)果與遞歸計算法中采用8 階Runge-Kutta 方法(ODE89)求解結(jié)果的對比,此時兩種方法的最終計算結(jié)果基本一致.

圖4 耦合計算法與遞歸計算法對比Fig.4 Comparison of RCM and CCM

圖5 表示在不同采樣頻率f s下分別由遞歸計算法和耦合計算法得到的范德波爾系統(tǒng)前7 階非線性輸出頻率響應(yīng)函數(shù)的計算結(jié)果,數(shù)值計算步長δt=1/f s.圖5 表明由遞歸計算法計算出的高階非線性輸出頻率響應(yīng)函數(shù)的結(jié)果受數(shù)值計算的采樣頻率影響較大,而由耦合計算法計算出的高階非線性輸出頻率響應(yīng)函數(shù)的結(jié)果受采樣頻率影響較小.這是因為四階Runge-Kutta 法(MATLAB,ODE45)是一種自適應(yīng)步長的數(shù)值計算方法,所以耦合計算法求解過程中的計算步長會在多個方程中自適應(yīng)調(diào)節(jié),使得累計誤差較小,最終計算結(jié)果受采樣頻率影響不大.遞歸計算法在求解廣義伴隨線性方程時,步長在各階廣義伴隨線性方程中單獨地自適應(yīng)調(diào)節(jié),所以計算高階非線性輸出頻率響應(yīng)函數(shù)時,遞歸計算法產(chǎn)生的累計誤差較大,最終計算結(jié)果受采樣頻率影響較大.

圖5 RCM 和CCM 中采樣頻率對NOFRFs 的影響Fig.5 Effect of sampling frequency on NOFRFs in RCM and CCM

圖6 表示在不同采樣頻率f s下,由遞歸計算法計算出的第i階非線性輸出頻率響應(yīng)函數(shù)的均方誤差值(mean square error,MSE),其計算公式為

圖6 采樣頻率對RCM 計算非線性輸出頻率響應(yīng)函數(shù)均方誤差的影響Fig.6 Effect of sampling frequency on the mean square error of NOFRFs calculated by RCM

圖6 的計算結(jié)果表明:隨著采樣頻率的增加,由遞歸計算法計算的高階非線性輸出頻率響應(yīng)函數(shù)的均方誤差先減小后增大,這是因為初始采樣頻率增加時,計算步長減小,每步計算產(chǎn)生的誤差減小,最終使計算結(jié)果的均方誤差減小,但繼續(xù)增加采樣頻率后,計算步數(shù)增加,計算結(jié)果的累計誤差反而增大.所以遞歸計算法的累計誤差不能直接通過減少數(shù)值計算的步長來降低,需要選擇合適的采樣頻率才能達(dá)到較好的計算效果.

為了進一步比較遞歸計算法和耦合計算法的差異,需考慮噪聲因素對兩種計算方法計算結(jié)果的影響.在遞歸計算法和耦合計算法求解出的各階非線性輸出響應(yīng)yi(t) 中加入一定量的白噪聲,圖7 表示在不同信噪比(signal noise ratio,SNR)的高斯白噪聲影響下,遞歸計算法(RCM)和耦合計算法(CCM)的均方誤差值變化曲線.噪聲影響下遞歸計算法的均方誤差值計算公式為

圖7 表明在相同信噪比條件下,耦合計算法計算結(jié)果的均方誤差值比遞歸計算法計算結(jié)果的均方誤差值更小;同時,隨著信噪比的增大,由遞歸計算法計算出的高階非線性輸出頻率響應(yīng)函數(shù)的均方根值不再遞減,而耦合計算法計算出的高階非線性輸出頻率響應(yīng)函數(shù)的均方根值隨著信噪比的增大仍在減小.所以相比遞歸計算法,耦合計算法在噪聲影響下的計算誤差更小.

此外,圖8(a)和圖8(b)分別表示用廣義伴隨線性方程的耦合計算法和最小二乘法計算出的范德波爾振子非線性參數(shù)CE對系統(tǒng)前7 階非線性輸出頻率響應(yīng)函數(shù)的影響,非線性參數(shù)CEi={0.5,1,1.5}CE.圖8 的結(jié)果表明:

圖8 非線性參數(shù)對NOFRFs 的影響Fig.8 Effect of nonlinear parameters on NOFRFs

(1) 非線性輸出頻率響應(yīng)函數(shù)可以表示非線性系統(tǒng)的頻域特征;

(2) 高階的非線性輸出頻率響應(yīng)函數(shù)可以表示系統(tǒng)非線性參數(shù)的變化;

(3) 高階的非線性輸出頻率響應(yīng)函數(shù)能量逐階減小.

圖8(a)表明,非線性參數(shù)對系統(tǒng)的1 階非線性輸出頻率響應(yīng)函數(shù)沒有影響,這是因為1 階的非線性輸出頻率響應(yīng)函數(shù)本質(zhì)上是非線性系統(tǒng)的線性派生系統(tǒng)頻率響應(yīng)函數(shù),它不受系統(tǒng)非線性參數(shù)的影響.而非線性參數(shù)對其他高階的非線性輸出頻率響應(yīng)函數(shù)的影響較大,主要表現(xiàn)為非線性參數(shù)變化時,高階非線性輸出頻率響應(yīng)函數(shù)比低階非線性輸出頻率響應(yīng)函數(shù)變化明顯,這表明高階的非線性輸出頻率響應(yīng)函數(shù)對非線性參數(shù)變化敏感程度更高,所以可利用高階非線性輸出頻率響應(yīng)函數(shù)對非線性參數(shù)變化的敏感性用于結(jié)構(gòu)損傷檢測及故障診斷.

圖8(b)表明最小二乘法無法準(zhǔn)確計算系統(tǒng)高階的非線性輸出頻率響應(yīng)函數(shù).當(dāng)非線性參數(shù)變化較小時,由最小二乘法計算出的非線性輸出頻率響應(yīng)函數(shù)已無法表征系統(tǒng)的非線性參數(shù)變化,而廣義伴隨線性方程的耦合計算方法比最小二乘方法計算高階非線性輸出頻率響應(yīng)函數(shù)的誤差小,最終計算結(jié)果更準(zhǔn)確.

2.2 達(dá)芬(Duffing)系統(tǒng)

以帶有立方剛度的非線性彈簧構(gòu)成的達(dá)芬振子為例,系統(tǒng)的NDE 模型為

式中,m為振子質(zhì)量,c為線性阻尼,k表示線性剛度,k3表示非線性剛度.

將式(44)代入式(26)中得系統(tǒng)的前五階廣義伴隨線性方程為

許多工程系統(tǒng)可以看作是單輸入單輸出的非線性系統(tǒng),其NDE 模型可用二階微分方程表示為[25]

式中m,c,k,kj,j=2,3,···,J為系統(tǒng)參數(shù),J為系統(tǒng)輸出非線性項的最高階數(shù),將式(5)代入式(46)得

Feijoo 等[25]指出式(47)的伴隨線性方程為

式中,j1+j2+···+jl=n.

將式(44) 代入式(48) 后可知達(dá)芬系統(tǒng)的前5 階伴隨線性方程與其廣義伴隨線性方程(45)的表達(dá)式完全一致,所以廣義伴隨線性方程法分析式(44)這類僅有輸出整數(shù)次冪非線性項的系統(tǒng)時,它與伴隨線性方程法的結(jié)果相同.同時,由伴隨線性方程式(47)和式(48)可知,伴隨線性方程法不能用于分析含有輸入輸出及其導(dǎo)數(shù)交叉項的非線性系統(tǒng)(例如范德波爾系統(tǒng)),所以廣義伴隨線性方程法比伴隨線性方程法的適用范圍更廣.

在零初始條件下,式(44)所表示的達(dá)芬系統(tǒng)第2 階與第4 階非線性輸出響應(yīng)y2(t),y4(t) 均為零解,即y2(t)=y4(t) .式(45)為線性微分方程,因為諧波函數(shù)是微分算子的本征函數(shù),所以對線性系統(tǒng)而言,輸入激勵頻率與系統(tǒng)響應(yīng)頻率保持一致,如u(t)=cos(2πft),系統(tǒng)第1 階非線性輸出響應(yīng)y1(t) 的頻率成分與系統(tǒng)輸入激勵u(t) 的頻率成分保持一致,即為f.第3 階非線性輸出y3(t) 的頻率成分與系統(tǒng)3 階廣義輸入激勵-k3y13(t) 的頻率成分相同,即為f,3f.同理,系統(tǒng)第5 階非線性輸出響應(yīng)y5(t) 的頻率成分與系統(tǒng)第5 階的廣義輸入激勵-3k3y12(t)y3(t) 的頻率成分相同,即為f,3f,5f.由于達(dá)芬系統(tǒng)的偶數(shù)階非線性輸出響應(yīng)y2k(t)=0,所以在諧波激勵下的達(dá)芬系統(tǒng)只產(chǎn)生了奇數(shù)次諧波分量,這是因為達(dá)芬系統(tǒng)式(44) 的非線性項為y3(t),它僅包含輸出y(t) 的3 次非線性項,即在式(26)中,只有輸出非線性項系數(shù)Cp,0(l1,l2,···,lp)≠0,其余高階廣義伴隨線性方程的非線性項系數(shù)均為0.由式(45)表示的達(dá)芬系統(tǒng)廣義伴隨線性方程也揭示了非線性系統(tǒng)在諧波激勵時“倍數(shù)頻率”現(xiàn)象產(chǎn)生的原因.

由式(45)得,非線性系統(tǒng)的廣義伴隨線性方程中線性算子 TL為

達(dá)芬系統(tǒng)的派生系統(tǒng)為

對單輸入單輸出的非線性系統(tǒng),廣義伴隨線性方程的線性算子與非線性系統(tǒng)的派生系統(tǒng)表達(dá)式相同.派生系統(tǒng)的無阻尼固有頻率.當(dāng)系統(tǒng)(44)輸入為諧波激勵u(t)=cos(2πfnt/3) 或激勵頻率f=fn/r,r=3,5,7,··· 時,由式(45)可知,達(dá)芬系統(tǒng)第3 階廣義伴隨線性方程的廣義輸入激勵-k3y13(t) 中產(chǎn)生了與派生系統(tǒng)固有頻率fn相同的諧波分量,此時系統(tǒng)第3 階非線性輸出響應(yīng)y3(t) 出現(xiàn)了“次諧共振”現(xiàn)象.

表3 表示達(dá)芬系統(tǒng)前5 階廣義伴隨線性方程的頻率分布,其中f為單諧波激勵的輸入頻率.

表3 廣義伴隨線性方程頻率分布Table 3 Frequency distribution of GALEs

設(shè)達(dá)芬系統(tǒng)式(44) 的參數(shù)為:m=1 kg,c=30 N·s/m,k=1.0×104N/m,k3=5.0×108N/m3,輸入u(t) 為

達(dá)芬振子線性派生系統(tǒng)的無阻尼固有頻率為穩(wěn)態(tài)響應(yīng)由Runge-Kutta 方法計算得到.圖9(a)～圖9(c)表示由廣義伴隨線性方程計算的前5 階非線性輸出響應(yīng)頻譜,其中圖9(a)表明第1 階非線性輸出響應(yīng)y1(t) 的頻率只有單倍頻成分 5 Hz .圖9(b)表明第3 階非線性輸出響應(yīng)y3(t) 的頻率包括單倍頻成分 5 Hz 和3 倍頻成分 15 Hz,圖9(c)表明第5 階非線性輸出響應(yīng)y5(t) 的頻率包括單倍頻成分5 Hz,3 倍頻成分 15 Hz 以及5 倍頻成分 25 Hz .由于諧波激勵u(t) 的頻率 5 Hz 接近達(dá)芬系統(tǒng)線性派生系統(tǒng)無阻尼固有頻率 16 Hz 的1/3,所以系統(tǒng)的第3 階與第5 階非線性輸出響應(yīng)在 15 Hz 處出現(xiàn)共振.

圖9 前5 階非線性輸出響應(yīng)頻譜Fig.9 Frequency spectrums of the first 5 orders of nonlinear output response

3 結(jié)論

本文基于NDE 模型表示的非線性系統(tǒng)廣義頻率響應(yīng)函數(shù)遞推公式及Volterra 級數(shù)理論,推導(dǎo)出了在一般激勵下由NDE 模型表示的非線性系統(tǒng)廣義伴隨線性方程計算公式.該公式表明非線性系統(tǒng)的輸出響應(yīng)可以用一系列線性方程組的解來表示,該結(jié)果的正確性通過仿真實驗研究進行了驗證.針對NDE 模型的廣義伴隨線性方程的數(shù)值求解問題,提出了耦合計算和遞歸計算兩種方法,研究結(jié)果表明,耦合計算法計算精度高,最終計算結(jié)果受數(shù)值計算采樣頻率的影響很小,能夠準(zhǔn)確地計算出系統(tǒng)非線性輸出頻率響應(yīng)函數(shù).此外,基于NDE 模型的廣義伴隨線性方程,研究了非線性參數(shù)對非線性輸出頻率響應(yīng)函數(shù)的影響,為基于非線性輸出頻率響應(yīng)函數(shù)在工程系統(tǒng)和結(jié)構(gòu)的狀態(tài)檢測及故障診斷中的應(yīng)用提供了一種新的思路.最后,以達(dá)芬振子系統(tǒng)為例,根據(jù)線性算子理論和廣義伴隨線性方程法分析了非線性系統(tǒng)中幾種典型非線性現(xiàn)象產(chǎn)生的理論依據(jù).本文研究表明,NDE 模型廣義伴隨線性方程法拓寬了伴隨線性方程法的適用范圍,對非線性系統(tǒng)的分析與設(shè)計具有重要意義.此外,由于NDE 模型廣義伴隨線性方程遞歸表達(dá)式計算復(fù)雜,不利于計算機編程實現(xiàn)及實際工程應(yīng)用,如何改進其計算方法需進一步研究.