張 波 張文博 彭志科 ,?,

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

? (上海交通大學(xué)機(jī)械系統(tǒng)與振動(dòng)國家重點(diǎn)實(shí)驗(yàn)室,上海 200240)

引言

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

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

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

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

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

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

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

微分算子 D 定義為

針對(duì)由式(1)所描述的系統(tǒng),當(dāng)系統(tǒng)輸入是具有實(shí)際物理意義的信號(hào)u(t) 時(shí),該系統(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)在零平衡點(diǎn)漸近穩(wěn)定,則該系統(tǒng)的輸出y(t) 可以用Volterra 級(jí)數(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 級(jí)數(shù)的截?cái)嚯A數(shù).對(duì)式(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ù)諧波探測(cè)法推導(dǎo)出了由NDE 模型表示的計(jì)算非線性系統(tǒng)廣義頻率響應(yīng)函數(shù)遞歸表達(dá)式

式中,LN=(l1,l2,···,lp+q) 為與Cp,q(l1,l2,···,lp+q) 相對(duì)應(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 階以后很難直觀表示出來,這對(duì)廣義頻率響應(yīng)函數(shù)在工程中的分析估計(jì)造成了極大困難.為了克服這個(gè)困難,Billings 等[7]提出了非線性系統(tǒng)的非線性輸出頻率響應(yīng)函數(shù)的概念.對(duì)于單輸入單輸出非線性系統(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 時(shí),Gn(jω)=G1(jω)表示線性系統(tǒng)的頻率響應(yīng)函數(shù).

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

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

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

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

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

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

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

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

式中

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

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

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

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

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

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

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

式中,n∈N?,LN=(l1,l2,···,lp+q) 是式(8)中與Cp,q(l1,l2,···,lp+q) 相對(duì)應(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) 與在時(shí)域中的類比,它揭示了一種關(guān)于非線性系統(tǒng)時(shí)域與頻域的新關(guān)系.

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

將式(27)代入式(29)可得Vn(t) 的表達(dá)式,其中廣義輸入Vn(t) 是由系統(tǒng)輸入激勵(lì)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 廣義伴隨線性方程計(jì)算非線性輸出頻率響應(yīng)函數(shù)的步驟

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

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

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

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

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

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

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

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

表2 耦合計(jì)算法Table 2 Coupled computational method

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

本節(jié)以范德波爾振子為例,利用數(shù)值仿真實(shí)驗(yàn)對(duì)上述理論分析結(jié)果的正確性進(jìn)行驗(yàn)證.針對(duì)廣義伴隨線性方程求解非線性輸出頻率響應(yīng)函數(shù)的數(shù)值計(jì)算問題,分析了遞歸計(jì)算法與耦合計(jì)算法的計(jì)算結(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)的輸入激勵(lì),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,輸入激勵(lì)

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

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

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

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

其中,αm∈R+,表示系統(tǒng)輸入激勵(lì)的次數(shù),N表示系統(tǒng)的截?cái)嚯A數(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 曲線為最小二乘法的計(jì)算結(jié)果.

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

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

圖4 耦合計(jì)算法與遞歸計(jì)算法對(duì)比Fig.4 Comparison of RCM and CCM

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

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

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

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

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

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

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

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

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

圖8(b)表明最小二乘法無法準(zhǔn)確計(jì)算系統(tǒng)高階的非線性輸出頻率響應(yīng)函數(shù).當(dāng)非線性參數(shù)變化較小時(shí),由最小二乘法計(jì)算出的非線性輸出頻率響應(yīng)函數(shù)已無法表征系統(tǒng)的非線性參數(shù)變化,而廣義伴隨線性方程的耦合計(jì)算方法比最小二乘方法計(jì)算高階非線性輸出頻率響應(yīng)函數(shù)的誤差小,最終計(jì)算結(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)輸出非線性項(xiàng)的最高階數(shù),將式(5)代入式(46)得

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

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

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

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

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

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

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

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

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