錢(qián)征文，程 禮，陳 衛(wèi)，李應(yīng)紅

(空軍工程大學(xué) 工程學(xué)院，西安 710038)

非線性振動(dòng)分析中一個(gè)重要方面是研究系統(tǒng)的周期解及其穩(wěn)定性。一般說(shuō)來(lái)，借助于數(shù)值方法可獲得系統(tǒng)的穩(wěn)定周期解，但由于其解曲線敏感依賴(lài)于初值的選擇，將很難獲得多解及其不穩(wěn)定解，因而需要借助于各種近似方法。

對(duì)于簡(jiǎn)單的具有多重解的振動(dòng)方程，如Duffing振子，許多作者［1－7]已利用諧波平衡法分析了其穩(wěn)定響應(yīng)，即假定振動(dòng)響應(yīng)可被表示成簡(jiǎn)單的截?cái)喔道锶~級(jí)數(shù)，從而將Duffing振子方程化為一個(gè)非線性代數(shù)方程組。但對(duì)于一些復(fù)雜的具有多重解的振動(dòng)方程，其對(duì)應(yīng)的非線性代數(shù)方程組非常復(fù)雜，可能依賴(lài)于不同的初值存在多重解，并且不穩(wěn)定解支的收斂域往往都很小，很難去直接求解。常用的求解非線性代數(shù)方程組的方法，如:牛頓迭代法、梯度下降法、擬牛頓法等對(duì)初值的選取是敏感依賴(lài)的，在求解多解問(wèn)題時(shí)效果不能令人滿意。此外，當(dāng)追蹤參數(shù)變化下系統(tǒng)的振動(dòng)響應(yīng)時(shí)，傳統(tǒng)的方法是在給定的初值下重復(fù)迭代求解過(guò)程，效率較低，初值選取不當(dāng)不僅浪費(fèi)很多計(jì)算時(shí)間，甚至有可能造成求解失敗。

針對(duì)上述問(wèn)題，考慮到同倫算法可以擴(kuò)大迭代格式的收斂域，將其引入到方程的求解中。同時(shí)為避免重復(fù)迭代求解效率低的問(wèn)題，采用預(yù)測(cè)校正算法來(lái)追蹤解曲線，從而確定非線性方程的穩(wěn)定周期解和不穩(wěn)定周期解。

1 求解方法推導(dǎo)

對(duì)于非線性振動(dòng)方程，一般表示為:

由諧波平衡法可以得到關(guān)于n個(gè)諧波系數(shù)以及外激勵(lì)頻率ω的非線性方程組:

式中，u=(a1，a2，…，an)，ai為第 i個(gè)諧波系數(shù)，n 的取值與選取的近似解的階數(shù)有關(guān)。

1.1 預(yù)測(cè)校正算法［8]

將非線性方程組(2)的解曲線問(wèn)題轉(zhuǎn)化為常微分方程的柯西問(wèn)題:

這樣，采用歐拉法或者其他更高階數(shù)值方法求得一條通過(guò)(u0，ω0)的積分曲線就是方程組(2)的解曲線，得到了諧波系數(shù)就可以確定非線性振動(dòng)方程(1)的近似解。

采用歐拉型積分公式，則有:

在式(3)的一步步求解過(guò)程中，可以不時(shí)地利用式(2)進(jìn)行一般牛頓迭代的修正:

使近似解點(diǎn)充分逼近解曲線，誤差滿足精度要求。

上述算法的具體實(shí)施過(guò)程中，對(duì)于有多重解存在的情況，步長(zhǎng)Δω的選取十分重要，較大的步長(zhǎng)會(huì)造成解突然從一個(gè)解支跳到另一個(gè)解支。為避免這種情況發(fā)生，在保證計(jì)算效率的同時(shí)可以盡量選取較小的步長(zhǎng)。

1.2 同倫算法

利用預(yù)測(cè)校正法求解非線性方程組(2)的過(guò)程中，需要確定初始值(u0，ω0)，通常情況下初始值的選取取決于計(jì)算者的經(jīng)驗(yàn)，并無(wú)特殊考慮，因此，初始值選取是方程求解的關(guān)鍵。而當(dāng)非線性方程組存在多解的情況時(shí)，初始值的選擇就更為困難。

同倫算法［9]是一種有效的方法，在擴(kuò)大迭代格式收斂域方面起著關(guān)鍵作用，理論上對(duì)迭代初始值的選取沒(méi)有任何限制。將同倫算法引入非線性方程組的求解中，就可以利用其全局搜索能力來(lái)提高解的收斂性。

為求非線性方程F(u，ω)=0的解，構(gòu)造一個(gè)帶有參數(shù) t的映射 H∶［0，1] × Rn→Rn，使得:

其中，方程 G(u，ω)=0的解已知，這個(gè)映射H［t，(u，ω)] 就稱(chēng)為同倫映射。

通過(guò)同倫映射，把求解方程F(u，ω)=0轉(zhuǎn)化為求解同倫方程H(t，(u，ω))的解。如果同倫方程對(duì)任何0≤t≤1有解u(t)存在，那么對(duì)應(yīng)于Rn中的解曲線u(t)，起點(diǎn)為u(0)，終點(diǎn)為u(1)，且有:

由式(7)可以看出，終點(diǎn)u(1)就是所求的方程F(u，ω)=0的解。

由于方程G(u，ω)=0的解是已知的，也就是解的初始點(diǎn)已知，追蹤解曲線的過(guò)程可以運(yùn)用前面的預(yù)測(cè)校正算法，從初始點(diǎn)起追蹤曲線u(t)直到終點(diǎn)u(1)，就可以得到原方程F(u，ω)=0的解。

常見(jiàn)的同倫映射有以下三種形式:

(1)不動(dòng)點(diǎn)同倫映射:H(t，x)=tF(x)+(1－t)·［x－x(0)] ;

(2)牛頓同倫映射:H(t，x)=F(x)－(1－t)F·［x(0)] ;

(3)線性同倫映射:H(t，x)=tF(x)+(1－t)·G(x)。

本文采用牛頓同倫映射算法。

1.3 計(jì)算步驟

將同倫算法與預(yù)測(cè)校正算法相結(jié)合，求解外激勵(lì)參數(shù)變化下非線性振動(dòng)方程(1)周期解的步驟如下:

① 對(duì)方程(1)用諧波平衡法得到關(guān)于n個(gè)諧波系數(shù)以及外激勵(lì)頻率ω的非線性方程組F(u，ω)=0;

② 構(gòu)造同倫映射，方程G(u)=0的解即是迭代初始值;

③ 在ω=ω0下，借助預(yù)測(cè)校正算法，利用式(4)和式(5)，以t作為外激勵(lì)參數(shù)，追蹤解曲線u(t)，終點(diǎn)u(1)就是所要求的方程F(u，ω0)=0的解;

④ ω=ω0+Δω，以u(píng)(t)為迭代初始值，以ω作為外激勵(lì)參數(shù)，借助預(yù)測(cè)校正算法，得到方程組F(u，ω)=0的解;

⑤ 改變?chǔ)?，重?fù)上一步，得到每一個(gè)ω下方程組F(u，ω)=0的解 u(ω);

⑥ 將u(ω)的帶入相應(yīng)的諧波系數(shù)得到方程(1)的周期解。

對(duì)于有多解存在的情況，在構(gòu)造同倫映射時(shí)，可以通過(guò)改變方程G(u)=0使迭代初始值在大范圍內(nèi)變化，從而搜索得到不同的解。

2 計(jì)算驗(yàn)證

一般地，工程結(jié)構(gòu)中梁的非線性振動(dòng)問(wèn)題以及許多電路模型均可表示為非線性Duffing振子的形式，其振動(dòng)方程為:

式中，c為阻尼系數(shù)，d為線性剛度系數(shù)，b為非線性剛度系數(shù)，f為外激勵(lì)載荷，ω為外激勵(lì)頻率。

取方程(1)的一階近似解為x=A cos(ωt+φ)，則其理論解［10]為:

取系統(tǒng)參數(shù)如下:c=0.28，d=1，b=2.5，f=0.82。由式(9)可得非線性Duffing振子(式(8))的幅頻曲線如圖1所示。

運(yùn)用本文的算法對(duì)方程(8)進(jìn)行求解，假設(shè)其一階近似解為 x=A1cos(ωt)+A2sin(ωt)，幅值利用諧波平衡法可得:

圖1 非線性Duffing振子的理論解幅頻曲線Fig.1 Amplitude-frequency curve of nonlinear Duffing oscillator

通過(guò)初值搜索可以得到非線性Duffing振子的三個(gè)解支，如圖2所示。

三個(gè)解支中，虛線所表示的解支是不穩(wěn)定的周期解，兩外兩個(gè)解支則為穩(wěn)定的周期解。從圖中可以看到，采用本文的方法計(jì)算得到的結(jié)果與理論近似解(圖1)是吻合的。

將四階龍格－庫(kù)塔法求得的周期解與使用本文方法求解的結(jié)果進(jìn)行對(duì)比，如圖3所示。系統(tǒng)參數(shù)取值為:c=0.28，d=1，b=2.5，f=0.82，ω =2.0。

從圖中可以看出，對(duì)于幅值較小的周期解，兩種方法的結(jié)果非常接近［圖3(a)] ;對(duì)于幅值較大的周期解，兩種方法的結(jié)果有一定的差距［圖3(b)] ，這主要是因?yàn)橐浑A近似解中忽略了高階信號(hào)，存在階段誤差。實(shí)際上，使用本文方法求解時(shí)，近似解中考慮的階數(shù)越高，兩者的誤差就會(huì)越小。圖3(c)是用本文方法求得的不穩(wěn)定周期解的相圖，而龍格－庫(kù)塔法無(wú)法求得不穩(wěn)定周期解，只能得到穩(wěn)定的周期解。

3 算例

拉桿轉(zhuǎn)子是各類(lèi)航空發(fā)動(dòng)機(jī)及大型燃?xì)廨啓C(jī)中常用的一種轉(zhuǎn)子結(jié)構(gòu)形式，由于盤(pán)與盤(pán)之間靠拉桿螺栓連接，并不是一個(gè)連續(xù)的整體，其轉(zhuǎn)子振動(dòng)特性非常復(fù)雜。本文通過(guò)合理地簡(jiǎn)化，構(gòu)建轉(zhuǎn)子的運(yùn)動(dòng)方程，利用上述方法對(duì)其振動(dòng)特性進(jìn)行研究。簡(jiǎn)化后轉(zhuǎn)子的力學(xué)模型如圖4所示。

圖4 簡(jiǎn)化后轉(zhuǎn)子力學(xué)模型Fig.4 Sketch of simplified mechanical model

盤(pán)簡(jiǎn)化為質(zhì)量分別是m1、m2的兩個(gè)剛性圓盤(pán)，偏心量分別為e1、e2，兩盤(pán)質(zhì)量偏心矢量夾角為φ。拉桿結(jié)構(gòu)簡(jiǎn)化為等效的抗彎彈簧，不計(jì)彈簧的質(zhì)量，具有非線性的剛度特征，其彈性回復(fù)力可表示為p(x)=k3x+。兩端分別與不計(jì)質(zhì)量的彈性軸相連，軸的彎曲剛度分別為k1、k2。軸、抗彎彈簧的阻尼系數(shù)分別為c1、c2、c3。

只考慮該轉(zhuǎn)子模型在x－y平面內(nèi)的彎曲振動(dòng)，忽略圓盤(pán)重力的作用，則圓盤(pán)的質(zhì)心運(yùn)動(dòng)方程可表示為:

由于轉(zhuǎn)子在x和y方向的剛度、阻尼相同，不考慮重力作用時(shí)兩個(gè)方向的振動(dòng)相同，本文只考慮x方向的振動(dòng)。

設(shè)方程(11)的穩(wěn)態(tài)近似解為:x1=A1sin(ωt)+A2cos(ωt)，x2=A3sin(ωt)+A4cos(ωt)，y1=A5sin(ωt)+A6cos(ωt)，y2=A7sin(ωt)+A8cos(ωt)。

其中，A1、A2、A3、A4、A5、A6、A7、A8為待定系數(shù)。將近似解代入方程(11)中，平衡正弦和余弦函數(shù)系數(shù)后可得:

其中:

圖5 不同參數(shù)下的振幅隨轉(zhuǎn)速變化曲線Fig.5 Vibration curve with different parameters

分別取不同的系統(tǒng)參數(shù)，對(duì)拉桿轉(zhuǎn)子的雙穩(wěn)態(tài)特性進(jìn)行仿真計(jì)算和分析，計(jì)算中，參數(shù)取值如下:

圖5是轉(zhuǎn)子x方向振動(dòng)幅值隨轉(zhuǎn)子轉(zhuǎn)速變化的曲線圖。紅線和藍(lán)線分別是盤(pán)1和盤(pán)2質(zhì)心的振幅曲線，橫坐標(biāo)表示的是相對(duì)轉(zhuǎn)速，100%時(shí)轉(zhuǎn)速為2 000 rad/s。

通過(guò)計(jì)算可知，不穩(wěn)定解的收斂域很小，常用的求解微分方程的方法無(wú)法得到不穩(wěn)定解。利用本文提出的求解方法，可以快速地求出兩個(gè)穩(wěn)態(tài)解(圖5中實(shí)線所示);對(duì)于不穩(wěn)定的解(圖5中虛線所示)，通過(guò)同倫算法擴(kuò)大收斂域后，可以準(zhǔn)確追蹤不穩(wěn)定的解曲線。仿真計(jì)算驗(yàn)證了該方法的有效性。

4 結(jié)論

本文提出的方法可以解決迭代初值難以確定的問(wèn)題，有效地求解非線性振動(dòng)方程的多重解。由于迭代初值的選取沒(méi)有限制，通過(guò)全局搜索，不但可以計(jì)算穩(wěn)定的周期解，而且不穩(wěn)定的周期解也可以求出。通過(guò)使用預(yù)測(cè)－校正算法，可以快速求解外激勵(lì)參數(shù)變化下的周期解。此外，由于該方法減少了迭代求解的計(jì)算量，在使用諧波平衡法時(shí)，可以選取任意階數(shù)的諧波，因此該方法可以求解任意精度的周期解。本文的研究成果可進(jìn)一步推廣于求解其他類(lèi)型的非線性方程。

［1] Ickens R E M.Comments on the method of harmonic balance［J] .Journal of Sound and Vibration，1984，94:456 －460.

［2] Hamdan M N，Burton T D.On the steady state response and stability of non-linear oscillators using harmonic balance［J] .Journal of Sound and Vibration，1993，166(2):255 －266.

［3] 李鵬松，吳柏生.達(dá)芬－諧波振子的改進(jìn)解析逼近解［J] .振動(dòng)與沖擊，2004，23(3):113 － 116.

［4] 李銀山，張善元，董青田，等.用兩項(xiàng)諧波法求解強(qiáng)非線性Duffing方程［J] .太原理工大學(xué)學(xué)報(bào)，2005，36(6):690－693.

［5] 胡 輝，郭源君，鄭敏毅.一個(gè)非線性奇異振子的諧波平衡解［J] .振動(dòng)與沖擊，2009，28(2):121 － 123.

［6] 高 陽(yáng)，王三民，劉曉寧.一種改進(jìn)的增量諧波平衡法及其在非線性振動(dòng)中的應(yīng)用［J] .機(jī)械科學(xué)與技術(shù)，2005，24(6):663－665.

［7] 楊紹普，申永軍，劉獻(xiàn)棟.基于增量諧波平衡法的齒輪系統(tǒng)非線性動(dòng)力學(xué)［J] .振動(dòng)與沖擊，2005，23(4):40 －42.

［8] 虞 烈，劉 恒.軸承－轉(zhuǎn)子系統(tǒng)動(dòng)力學(xué)［M] .西安:西安交通大學(xué)出版社，2000.

［9] 何哲明，李 贊.同倫算法的常微分方程數(shù)值解法及在機(jī)構(gòu)學(xué)中的應(yīng)用［J] .機(jī)床與液壓，2003，1:214 －216.

［10] 郝亦清，李翠英.非線性振動(dòng)分析［M] .北京:北京理工大學(xué)出版社，1996.