非線性振動方程多重解求解方法

錢征文，程 禮，陳 衛，李應紅

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

非線性振動分析中一個重要方面是研究系統的周期解及其穩定性。一般說來，借助于數值方法可獲得系統的穩定周期解，但由于其解曲線敏感依賴于初值的選擇，將很難獲得多解及其不穩定解，因而需要借助于各種近似方法。

對于簡單的具有多重解的振動方程，如Duffing振子，許多作者［1－7]已利用諧波平衡法分析了其穩定響應，即假定振動響應可被表示成簡單的截斷傅里葉級數，從而將Duffing振子方程化為一個非線性代數方程組。但對于一些復雜的具有多重解的振動方程，其對應的非線性代數方程組非常復雜，可能依賴于不同的初值存在多重解，并且不穩定解支的收斂域往往都很小，很難去直接求解。常用的求解非線性代數方程組的方法，如:牛頓迭代法、梯度下降法、擬牛頓法等對初值的選取是敏感依賴的，在求解多解問題時效果不能令人滿意。此外，當追蹤參數變化下系統的振動響應時，傳統的方法是在給定的初值下重復迭代求解過程，效率較低，初值選取不當不僅浪費很多計算時間，甚至有可能造成求解失敗。

針對上述問題，考慮到同倫算法可以擴大迭代格式的收斂域，將其引入到方程的求解中。同時為避免重復迭代求解效率低的問題，采用預測校正算法來追蹤解曲線，從而確定非線性方程的穩定周期解和不穩定周期解。

1 求解方法推導

對于非線性振動方程，一般表示為:

由諧波平衡法可以得到關于n個諧波系數以及外激勵頻率ω的非線性方程組:

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

1.1 預測校正算法［8]

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

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

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

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

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

上述算法的具體實施過程中，對于有多重解存在的情況，步長Δω的選取十分重要，較大的步長會造成解突然從一個解支跳到另一個解支。為避免這種情況發生，在保證計算效率的同時可以盡量選取較小的步長。

1.2 同倫算法

利用預測校正法求解非線性方程組(2)的過程中，需要確定初始值(u0，ω0)，通常情況下初始值的選取取決于計算者的經驗，并無特殊考慮，因此，初始值選取是方程求解的關鍵。而當非線性方程組存在多解的情況時，初始值的選擇就更為困難。

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

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

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

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

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

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

常見的同倫映射有以下三種形式:

(1)不動點同倫映射: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 計算步驟

將同倫算法與預測校正算法相結合，求解外激勵參數變化下非線性振動方程(1)周期解的步驟如下:

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

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

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

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

⑤ 改變ω，重復上一步，得到每一個ω下方程組F(u，ω)=0的解 u(ω);

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

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

2 計算驗證

一般地，工程結構中梁的非線性振動問題以及許多電路模型均可表示為非線性Duffing振子的形式，其振動方程為:

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

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

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

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

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

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

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

將四階龍格－庫塔法求得的周期解與使用本文方法求解的結果進行對比，如圖3所示。系統參數取值為:c=0.28，d=1，b=2.5，f=0.82，ω =2.0。

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

3 算例

拉桿轉子是各類航空發動機及大型燃氣輪機中常用的一種轉子結構形式，由于盤與盤之間靠拉桿螺栓連接，并不是一個連續的整體，其轉子振動特性非常復雜。本文通過合理地簡化，構建轉子的運動方程，利用上述方法對其振動特性進行研究。簡化后轉子的力學模型如圖4所示。

圖4 簡化后轉子力學模型Fig.4 Sketch of simplified mechanical model

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

只考慮該轉子模型在x－y平面內的彎曲振動，忽略圓盤重力的作用，則圓盤的質心運動方程可表示為:

由于轉子在x和y方向的剛度、阻尼相同，不考慮重力作用時兩個方向的振動相同，本文只考慮x方向的振動。

設方程(11)的穩態近似解為: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為待定系數。將近似解代入方程(11)中，平衡正弦和余弦函數系數后可得:

其中:

圖5 不同參數下的振幅隨轉速變化曲線Fig.5 Vibration curve with different parameters

分別取不同的系統參數，對拉桿轉子的雙穩態特性進行仿真計算和分析，計算中，參數取值如下:

圖5是轉子x方向振動幅值隨轉子轉速變化的曲線圖。紅線和藍線分別是盤1和盤2質心的振幅曲線，橫坐標表示的是相對轉速，100%時轉速為2 000 rad/s。

通過計算可知，不穩定解的收斂域很小，常用的求解微分方程的方法無法得到不穩定解。利用本文提出的求解方法，可以快速地求出兩個穩態解(圖5中實線所示);對于不穩定的解(圖5中虛線所示)，通過同倫算法擴大收斂域后，可以準確追蹤不穩定的解曲線。仿真計算驗證了該方法的有效性。

4 結論

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

［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] 李鵬松，吳柏生.達芬－諧波振子的改進解析逼近解［J] .振動與沖擊，2004，23(3):113 － 116.

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

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

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

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

［8] 虞 烈，劉 恒.軸承－轉子系統動力學［M] .西安:西安交通大學出版社，2000.

［9] 何哲明，李 贊.同倫算法的常微分方程數值解法及在機構學中的應用［J] .機床與液壓，2003，1:214 －216.

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