黃建亮 張兵許 陳樹輝

(中山大學(xué)應(yīng)用力學(xué)與工程系,廣州 510275)

引言

增量諧波平衡法(incremental harmonic balance,IHB) 是一個(gè)半解析半數(shù)值的方法,由Lau 和Cheung[1]于1981 年提出以來,是把增量過程牛頓-拉弗森(Newton-Raphson)迭代過程和諧波平衡過程(伽遼金(Galerkin)平均過程)兩者有機(jī)地結(jié)合起來,已被成功地推廣應(yīng)用于非線性振動(dòng)的各個(gè)領(lǐng)域.Lau[2]和陳樹輝[3]分別于1995 年和2007 年對(duì)該方法的發(fā)展及其應(yīng)用做了很好的綜述.在文獻(xiàn)[3-4]中,IHB 法已被總結(jié)為具有兩大優(yōu)點(diǎn):(1)應(yīng)用于分析強(qiáng)非線性振動(dòng)問題,并得到高精度解;(2)是一類求解非線性微分方程的計(jì)算方法,適用于大范圍參數(shù)變化的定量分析.

IHB 法已廣泛應(yīng)用在工程中的各個(gè)領(lǐng)域,包含航空航天、車輛交通、海洋工程、機(jī)械等各種非線性振動(dòng)問題的分析中[5-8].Wang 和Zhu[9]建立了非圓鏈輪的汽車皮帶傳動(dòng)系統(tǒng)的動(dòng)力學(xué)模型,用IHB法研究減小凸輪軸角度變化的效果.Mitra 等[10]利用IHB 法分析了強(qiáng)非線性半潛式平臺(tái)系統(tǒng)的主諧波響應(yīng)和高階次諧波響應(yīng).在機(jī)械振動(dòng)中,齒輪的振動(dòng)被廣泛研究,Li 等[11]分析了具有動(dòng)態(tài)齒隙的齒輪副系統(tǒng)在內(nèi)外周期激勵(lì)下的非線性動(dòng)力學(xué)特性;Zhou 等[12]基于IHB 法分析了考慮齒隙、阻尼、傳動(dòng)誤差和嚙合剛度的直齒圓柱齒輪副動(dòng)力學(xué)模型,發(fā)現(xiàn)了跳躍不連續(xù)現(xiàn)象和多個(gè)穩(wěn)定解共存.IHB 法對(duì)各種動(dòng)力學(xué)問題具有良好的求解效果.Wang 等[13]將非線性波傳播問題的控制方程轉(zhuǎn)換為相應(yīng)的具有多個(gè)時(shí)滯的非線性時(shí)滯微分方程,并應(yīng)用IHB 法進(jìn)行計(jì)算,大大簡(jiǎn)化推導(dǎo)的過程和計(jì)算的時(shí)間.Cheng 等[14]應(yīng)用IHB 法分析了附加恒力對(duì)準(zhǔn)零剛度隔振器的主共振和1/3 次諧波共振的影響.Li 和Chen[15]利用IHB 法研究了空氣壓縮機(jī)模型的周期運(yùn)動(dòng),混沌運(yùn)動(dòng)和分岔行為.Karli?i?等[16]利用IHB 法分析了基于馬蒂厄-達(dá)芬(Mathieu-Duffing)非線性振子模型的非線性壓電能量采集器的穩(wěn)態(tài)響應(yīng).Shen 等[17]用IHB 法研究了馬蒂厄-達(dá)芬振子的分岔和混沌路徑,識(shí)別出了一系列振子的倍周期分岔點(diǎn),發(fā)現(xiàn)它們近似服從普適標(biāo)度律.IHB 法也被應(yīng)用于分?jǐn)?shù)階非線性方程的計(jì)算過程中,Shen 等[18]將IHB 法應(yīng)用到分?jǐn)?shù)階非線性振子的動(dòng)力學(xué)分析中,建立了該振子周期解的一般形式來獲得高精度的解.對(duì)于振動(dòng)系統(tǒng)里含有多個(gè)不可約頻率的準(zhǔn)周期運(yùn)動(dòng)問題,Huang 等[19-21]將IHB 法發(fā)展成為多時(shí)間尺度IHB 法,應(yīng)用于求解含邊頻帶的準(zhǔn)周期運(yùn)動(dòng),自動(dòng)跟蹤非線性系統(tǒng)的準(zhǔn)周期運(yùn)動(dòng)頻率響應(yīng).竇蘇廣等[22]和蘇鸞鳴和葉敏[23]將IHB 法應(yīng)用于非線性系統(tǒng)參數(shù)識(shí)別,分析了在周期、分岔和混沌等復(fù)雜運(yùn)動(dòng)狀態(tài)下的參數(shù)識(shí)別,并驗(yàn)證了IHB 法相對(duì)于其他方法有明顯的抗噪性能.此外,他們還將增量諧波平衡非線性識(shí)別理論運(yùn)用到實(shí)驗(yàn)建模方法中[24-25].

為了提高IHB 法適應(yīng)性,人們對(duì)它進(jìn)行了各種改進(jìn),主要集中在兩個(gè)方面.第一方面是要提高IHB 法的計(jì)算效率.Wang 和Zhu[26]對(duì)非線性代數(shù)方程的余項(xiàng)使用快速傅里葉變換逼近,非線性代數(shù)方程的雅可比矩陣使用Broyden 擬牛頓迭代方法逼近,從而提高了諧波平衡過程的收斂速度;第二方面是解決系統(tǒng)響應(yīng)曲線在峰值(或拐點(diǎn))處的自動(dòng)跟蹤問題.為了可使IHB 法能更好的追蹤響應(yīng)隨參數(shù)變化情況,特別是在峰值(或拐點(diǎn))處自動(dòng)跟蹤,避免出現(xiàn)迭代不收斂的問題,可利用弧長(zhǎng)延拓法對(duì)IHB 法進(jìn)行改進(jìn)[4,26-27].另外,Zheng 等[28]在線性增量過程中引入Tikhonov 正則化來解決迭代中的不適定問題,從而改變IHB 法的收斂性能.然而對(duì)IHB 法的各種改進(jìn),沒有涉及到如何選擇初值條件的問題,在給定的初值條件有時(shí)并不能保證IHB 法的收斂性.

在IHB 法的牛頓-拉弗森迭代過程中,初值可以任意設(shè)定,實(shí)際上IHB 法的牛頓-拉弗森迭代過程相當(dāng)于對(duì)初值進(jìn)行多次嘗試,并不能保證初值的收斂與否,當(dāng)?shù)恼`差增大時(shí)相當(dāng)于換一個(gè)初值進(jìn)行重新迭代嘗試.本文針對(duì)這一問題,引入回溯線搜索(backtracking line search,BLS)優(yōu)化算法[29]來調(diào)節(jié)IHB 法的收斂步長(zhǎng),使步長(zhǎng)逐漸減小直至滿足弱沃爾夫(Wolfe)條件,即下降性條件,來對(duì)初值進(jìn)行迭代計(jì)算.在此基礎(chǔ)上,引入迭代負(fù)梯度方向,再利用迭代優(yōu)化算法狗腿法(dogleg method)[30]的思想,加快其收斂速度,提高了初值嘗試的局部收斂范圍,從而提高了IHB 法牛頓-拉弗森迭代的收斂性.

1 傳統(tǒng)的增量諧波平衡法

考慮一個(gè)一般形式的n維非線性振動(dòng)系統(tǒng),可表示為如下平衡微分方程

其中,“·”表示對(duì)時(shí)間t的微分,xi(t) 是位移變量,pi表示第i個(gè)外激勵(lì),Φi是二階的非線性常微分方程.引入新的時(shí)間變量 τ,定義為

其中,ω 為系統(tǒng)的振動(dòng)響應(yīng)頻率.方程(1)可表示為

其中,“ ′ ”表示對(duì)時(shí)間 τ 的微分.

IHB 法的第一步是增量過程,即牛頓-拉弗森迭代過程.設(shè)xi0,pi0和 ω0是方程(3)的解,則其鄰近解可用增量形式表示為

其中,Δxi,Δpi和Δ ω 為增量.把方程(4) 代入方程(3),并略去 Δxi,Δpi和 Δ ω 的高階小量,得到了線性化的增量方程

其中,下標(biāo)0 表示為函數(shù)增量的初值.方程(5)可表示為矩陣形式

IHB 法的第二步是諧波平衡過程,即伽遼金平均過程.將n維非線性系統(tǒng)的穩(wěn)態(tài)解及其增量都展開成傅里葉級(jí)數(shù)

其中,nc和ns是正整數(shù),表示諧波項(xiàng)項(xiàng)數(shù);Cs,Ai和ΔAi分別表示為

那么向量X0及其增量可分別表示為傅里葉系數(shù)向量A及其增量ΔA

積分上式并整理歸并為以 ΔA,Δ ω 和 ΔP為未知量的代數(shù)方程組

這里,KA和Rp是nt×nt的矩陣,R和Rω是nt×1 的列向量,nt=n×(nc+ns) .當(dāng)要分析某一固定外激勵(lì)幅值下的頻率響應(yīng)曲線時(shí),在方程(14)中取P為固定值,即 ΔP=0,那么方程(14)變成

當(dāng)要分析某一固定外激勵(lì)頻率下的各響應(yīng)幅值隨外激勵(lì)幅值變化的響應(yīng)曲線時(shí),在方程(14)中取ω為固定值,即 Δ ω=0,那么方程(14)變成

在用傳統(tǒng)的IHB 法對(duì)方程式(19)或式(20)的求解過程中,要給定一個(gè)假定的初值,然后可先利用諧波平衡過程再利用增量過程,或也可先利用增量過程再利用諧波平衡過程,追蹤出所有的解.在增量過程中,可應(yīng)用頻率增量,振幅增量,或弧長(zhǎng)增量,具體可參見文獻(xiàn)[3-4].

如上所述,在用傳統(tǒng)IHB 法求解過程中,事先要給定一假定的初值,因此涉及到給定初值的選取問題.一般說來,當(dāng)碰到這類問題時(shí),一般是用三種途徑來解決:(1)憑IHB 法使用者的經(jīng)驗(yàn),事先給定一接近于精確解作為IHB 法的初值;(2)以線性化系統(tǒng)得到的解作為IHB 法的初值;(3)利用攝動(dòng)法等求得小參數(shù)范圍內(nèi)的弱非線性振動(dòng)的解作為IHB 法的初值.在傳統(tǒng)IHB 法用牛頓-拉弗森迭代過程中,如給定的初值在通往收斂的路徑上,那么解就會(huì)收斂至精確解,即此初值點(diǎn)是在其收斂范圍內(nèi).為了擴(kuò)大初值點(diǎn)的選擇范圍,在用牛頓-拉弗森迭代過程中,本文將對(duì)迭代的增量進(jìn)行兩個(gè)方面的改進(jìn).

2 引入BLS 優(yōu)化算法的改進(jìn)增量諧波平衡法(GIHB1)

其中,向量KAk和Rk分別表示KA和R在第k次迭代的結(jié)果.

現(xiàn)對(duì)傳統(tǒng)IHB 法中的迭代增量進(jìn)行改進(jìn).首先,引入目標(biāo)函數(shù)

當(dāng) ‖R‖→0 時(shí),f(A)→0 .再引入一次函數(shù)和基于Cauchy-Point 的二次函數(shù)

其中,向量Ak表示A在第k次迭代的結(jié)果.

當(dāng)初值選擇恰當(dāng)時(shí),‖ ΔA‖ 越小,對(duì)f(Ak+ΔAk)在點(diǎn)Ak上進(jìn)行模擬,m2(Ak+ΔAk) 與f(Ak+ΔAk) 的差越小.在迭代過程中,用m2(Ak+ΔAk) 代替f(Ak+ΔAk)進(jìn)行求解,則增量變?yōu)?/p>

當(dāng)KAk為非奇異矩陣時(shí),等價(jià)于原牛頓-拉弗森迭代增量式(21),即使得迭代一步就到達(dá)m2(Ak+ΔAk)最小點(diǎn).

當(dāng)初值選擇不當(dāng)時(shí),可導(dǎo)致 ‖ ΔA‖ 較大,m2(Ak+ΔAk) 代替f(Ak+ΔAk) 吻合度較差,因此導(dǎo)致f(Ak+ΔAk) 的值不降反增,從而導(dǎo)致迭代不收斂.對(duì)此,引入BLS 優(yōu)化算法將其迭代步長(zhǎng)進(jìn)行改進(jìn),以為牛頓-拉弗森迭代方向,由于

因此牛頓-拉弗森迭代的方向 β1與負(fù)梯度方向呈銳角,滿足方向下降條件.在 β1前乘以系數(shù) α ,當(dāng)α較小時(shí),滿足弱Wolfe 條件時(shí)

即可選定迭代步長(zhǎng),然后再進(jìn)行IHB 法的第二步諧波平衡過程(Galerkin procedure),其中,常數(shù)c通常較小,一般可取c=1.0×10-4,該優(yōu)化算法的迭代過程如圖1 所示.在圖1 的算法中,為了保證較快的迭代速度,取 α =1 為牛頓-拉弗森迭代的初值,ρ在每次迭代時(shí)與 α 相乘,為了保證步α β1逐漸減小,取 0 ＜ρ ＜1 .隨著迭代次數(shù)的增加,α 值逐漸減小,即A的增量 α β1逐漸減小,直至滿足弱Wolfe 條件,然后跳出循環(huán),以得到Ak為初值,重新執(zhí)行牛頓-拉弗森迭代過程.

圖1 結(jié)合Backtracking line search (BLS)優(yōu)化算法的IHB 法(GIHB1)Fig.1 Algorithm for a generalized IHB method with backtracking line search (BLS) method (GIHB1)

3 引入狗腿法思想并結(jié)合BLS 算法的改進(jìn)IHB 法(GIHB2)

為了加快BLS 算法的迭代運(yùn)算速度,應(yīng)用狗腿法的思想,在此過程中加入兩個(gè)調(diào)節(jié)參數(shù)用于控制迭代方向.當(dāng) ‖ ΔA‖ 較小時(shí),為負(fù)梯度方向,即最快的下降方向.m2(Ak+ΔAk) 在 β2方向上的最小值對(duì)應(yīng)的系數(shù) αmin可由

系數(shù) αmin的主要作用是調(diào)節(jié)負(fù)梯度方向的步長(zhǎng),使得該步長(zhǎng)和模擬函數(shù)牛頓-拉弗森迭代的步長(zhǎng)處于同一個(gè)量級(jí).在狗腿法思想基礎(chǔ)上,引入兩個(gè)調(diào)節(jié)參數(shù) α1和 α2用于加快BLS 算法的迭代運(yùn)算速度.令Ψ1k=α1αminβ2,Ψ2k=β1,其中,α1是調(diào)節(jié)目標(biāo)函數(shù)f(Ak) 和二次模擬函數(shù)m2(Ak+ΔAk) 的迭代方向,α1∈(0,0.1) .當(dāng) ΔA沿 著 0 →Ψ1k→Ψ2k方向時(shí),m2(Ak+ΔAk) 逐 漸下降到最小值,其中 0 表示為零向量.

引入函數(shù)

用于調(diào)節(jié)迭代的方向,其中,α2起線性化分割的作用,其值可取 α2∈(0,0.1) .同時(shí),引入m2(Ak+ΔAk)與f(Ak+ΔAk) 的擬合系數(shù) ρk,可表示為

此系數(shù) ρk越大時(shí),m2(Ak+ΔAk) 與f(Ak+ΔAk)的吻合程序越高,f(Ak+ΔAk) 越滿足下降條件.算法的迭代過程可利用系數(shù) ρk值的大小來判斷 α 是否滿足條件.對(duì)于系數(shù) ρk界定一個(gè)閾值 ρ0,當(dāng) ρk＞ρ0時(shí),系數(shù) α 滿足條件;否則繼續(xù)減小,直至滿足ρk＞ρ0,該算法的迭代過程如圖2 所示.

在圖2 所示的算法中,α =1 時(shí)迭代方式為牛頓-拉弗森迭代,α 與F(α) 的模長(zhǎng)有正相關(guān)的關(guān)系.隨著迭代次數(shù)的增加,α 值逐漸減小,即A的增量F(α) 的模長(zhǎng)逐漸減小,直至 ρk＞ρ0,然后跳出循環(huán),得到初值A(chǔ)k,然后重新執(zhí)行牛頓-拉弗森迭代過程.

圖2 結(jié)合BLS 算法和狗腿算法的改進(jìn)IHB 法(GIHB2),其中在狗腿算法的函數(shù) F (α) 中引入2 個(gè)調(diào)節(jié)參數(shù) α 1 和α2Fig.2 Algorithm for a generalized IHB method (GIHB2) combined with BLS method and dogleg method that contains two parametersα1 and α2 in the functionF(α)

值得注意的是上述兩種改進(jìn)的IHB 法隨著 α 的減小,迭代效率會(huì)逐漸降低,當(dāng) α ＜10-6時(shí),即認(rèn)為迭代不收斂.為了提高收斂性的同時(shí),提高迭代的效率,在算法實(shí)現(xiàn)過程中,可設(shè) α 的下界值為0 .01,若α ＜0.01 ,則令A(yù)k+1=Ak+0.1β1改變初值,重新執(zhí)行算法的迭代過程.

4 算例

利用3 種IHB 法,即傳統(tǒng)的IHB 法、結(jié)合BLS 算法的GIHB1 法和引入狗腿思想并結(jié)合BLS 算法的GIHB2 法,對(duì)兩個(gè)算例的初值問題進(jìn)行分析,一個(gè)是經(jīng)典的單自由度范德波爾(van der Pol)振蕩,另一個(gè)是耦合的兩自由度范德波爾振蕩.在IHB 法的求解過程中,此二類范德波爾振蕩的極限環(huán)周期解依賴于給定的初值條件.為了檢驗(yàn)上述3 種IHB 法計(jì)算結(jié)果的準(zhǔn)確性,本文采用四階龍格-庫(kù)塔(Runge-Kutta,RK)數(shù)值法作為比較.四階R-K數(shù)值法是直接對(duì)微分方程采用時(shí)間積分法,它可以較準(zhǔn)確地給出系統(tǒng)某一時(shí)刻的位移、速度和加速度的數(shù)值.

4.1 van der Pol 振蕩

單自由度van der Pol 振蕩方程為

其中,“·”表示對(duì)間t的導(dǎo)數(shù),參數(shù) λ 在區(qū)間 [ 0,1.0] 內(nèi)取值用來分析系統(tǒng)響應(yīng)的霍普夫分岔.當(dāng)單自由度van der Pol 振蕩方程(32)中的2 個(gè)參數(shù) ε 和 λ 確定時(shí),該振蕩方程具有唯一一個(gè)極限環(huán).這類van der Pol 振蕩方程在航空航天、力學(xué)以及電子等工程領(lǐng)域都有很廣泛的應(yīng)用.例如,Motta 和Malzacher[31]發(fā)展了湍流經(jīng)過機(jī)翼面的van der Pol 模型,根據(jù)此模型構(gòu)建了開環(huán)與閉環(huán)的流量控制.Akhtar 等[32]提出了水動(dòng)力作用在圓柱/橢圓柱結(jié)構(gòu)上的范德波爾-達(dá)芬振蕩模型,用此模型得到的計(jì)算結(jié)果與CFD 的結(jié)果一致.Facchinetti 等[33]基于van der Pol 振蕩方程研究了在渦振中結(jié)構(gòu)與波的耦合振蕩.

值得注意的是,經(jīng)典攝動(dòng)法只能適用于 ε 為小參數(shù)時(shí)van der Pol 振蕩方程的求解,而不適用于大參數(shù)時(shí)的求解.Burton[34]說明了當(dāng) ε ＞1.0 時(shí),利用攝動(dòng)法去求解會(huì)得到錯(cuò)誤的結(jié)果.而IHB 法的優(yōu)勢(shì)是不僅適用于弱非線性系統(tǒng),也適用于強(qiáng)非線性系統(tǒng),不失一般性,文中 ε 可取大參數(shù)作為分析.

引入新的時(shí)間變量

方程(33)變?yōu)?/p>

其中,“ ′ ”表示為時(shí)間變量 τ 的導(dǎo)數(shù).

因方程(34)只含奇次項(xiàng),系統(tǒng)是對(duì)稱的振蕩響應(yīng),可將x展開為含有m個(gè)奇次諧波項(xiàng)的傅里葉級(jí)數(shù)

根據(jù)式(14),得

其中,KA是 2m×2m的矩陣,對(duì)應(yīng)于傅里葉系數(shù)向量

向量R和Rω可分別由式(16)和式(17)計(jì)算得到.在本文中取誤差向量R的模 ‖R‖ 小于1 .0×10-7作為IHB 法諧波平衡過程的收斂條件.

顯然,式(36)中所含的未知量的數(shù)目比方程的數(shù)目多2 個(gè),由于系統(tǒng)響應(yīng)頻率 ω 是未知的,因此必須選定 Δ λ 為主動(dòng)增量,再選定式(37)中的某一諧波項(xiàng)系數(shù)為參考量,如b1,那么該參考量的增量為零,即 Δb1=0,這樣方程(36)就可求解.具體步驟可參見文獻(xiàn)[3-37]等.

現(xiàn)利用上述3 種方法分別對(duì)van der Pol 振蕩方程(32)求解,其中,λ =0.87 ,ε =5.0 ,選定b1為參考量.同時(shí),在式(35)中取諧波項(xiàng)項(xiàng)數(shù)為m=25 .圖3所示為在取不同初值時(shí)用3 種方法得到的誤差‖R‖與迭代次數(shù)的關(guān)系.從圖中可以得到:

(1)圖3(a)中選定初值為頻率 ω0=0.58,諧波項(xiàng)系數(shù)a1=1.0 ,b1=1.0,其余諧波項(xiàng)系數(shù)設(shè)為零,可以看出傳統(tǒng)的IHB 無法收斂,另外兩種改進(jìn)方法均可收斂,GIHB1 法經(jīng)25 步迭代達(dá)到‖R‖＜1.0×10-7的收斂條件,GIHB2 法經(jīng)23 步迭代達(dá)到收斂條件,其響應(yīng)頻率為 ω =0.5877 ;

圖3 不同初值條件下用3 種IHB 法求解van der Pol 振蕩方程(32)迭代收斂情況,其中,ε =5.0 ,λ =0.87 .初值分別為:(a) ω 0=0.58 ,a 1=1.0 ,b 1=1.0 ;(b) ω 0=0.58,a 1=1.0,b 1=1.4 ;(c) ω 0=0.58,a 1=1.8,b1=-0.7Fig.3 Iteration convergence of solution of van der Pol Eq.(32) with different initial values using the three IHB methods,where ε =5.0,λ=0.87 .Initial values are (a) ω 0=0.58 ,a 1=1.0 ,b 1=1.0 ;(b) ω 0=0.58 ,a 1=1.0,b 1=1.4 ;(c) ω 0=0.58 ,a 1=1.8,b1=-0.7

(2)圖3(b)中選定初值為頻率 ω0=0.58,諧波項(xiàng)系數(shù)a1=1.0 ,b1=1.4,其余諧波項(xiàng)系數(shù)設(shè)為零,可以看出三種方法都可收斂,傳統(tǒng)IHB 法誤差 ‖R‖ 先增大再減少,其最大值可達(dá)到3113,然后開始下降,在誤差增大過程中,增量的模長(zhǎng)較大,易出現(xiàn)跳出牛頓-拉弗森迭代的現(xiàn)象,導(dǎo)致整個(gè)迭代不收斂.而對(duì)改進(jìn)后的兩種方法來說,誤差 ‖R‖ 可一直保持減少,直至滿足迭代條件.可以看出傳統(tǒng)IHB 法需要28 步迭代達(dá)到收斂條件,而改進(jìn)后的兩種GIHB 法分別只需要25 步和15 步迭代達(dá)到收斂條件;

(3)圖3(c)中選定初值為頻率 ω0=0.58,諧波項(xiàng)系數(shù)a1=1.8 ,b1=-0.7,其余諧波項(xiàng)系數(shù)設(shè)為零,可以看出傳統(tǒng)IHB 法無法收斂,兩種改進(jìn)GIHB 法均可收斂,GIHB1 法經(jīng)214 步迭代才能達(dá)到收斂條件,而GIHB2 法只需58 步迭代即可達(dá)到收斂條件.

圖4 所示為用3 種IHB 法求解van der Pol 振蕩方程(32) 的初值a1和b1的選定范圍,其他初值為ω0=0.58 及諧波項(xiàng)系數(shù)除值a1和b1外都為零.可以看出,兩種改進(jìn)后的方法其初值的選定范圍一樣,且均多于傳統(tǒng)IHB 法的初值選定范圍.圖5 所示為取λ=0.87,ε =5.0 時(shí)van der Pol 振蕩方程(32)的極限環(huán),可以看出3 種IHB 法求出的結(jié)果是一致的,是因?yàn)閮煞N改進(jìn)的GIHB 法只涉及到增量過程的改進(jìn),而諧波平衡過程沒有改變;同時(shí)3 種IHB 法求得的結(jié)果與用四階R-K 法求得的數(shù)值結(jié)果高度吻合.

圖4 用3 種IHB 法求得van der Pol 振蕩方程(32)解的初值 a1 和b1的選定范圍,其他初值為 ω 0=0.58 及其他諧波項(xiàng)系數(shù)設(shè)為零Fig.4 The regions for the initial value for a1 and b1 of solutions of van der Pol Eq.(32) using the three IHB methods,the other initial value for frequncy is ω 0=0.58 and the other initial harmonic terms are set to be zero

圖5 用3 種IHB 法和四階龍格-庫(kù)塔數(shù)值法求得范德波爾振蕩方程(32)的解,其中 λ =0.87 和ε=5.0Fig.5 Solutions of van der Pol Eq.(32) by the three IHB methods and numerical integration using the fourth-order R-K method,where λ=0.87 andε=5.0

4.2 耦合van der Pol 振蕩

眾多學(xué)者對(duì)工程中出現(xiàn)的耦合van der Pol 振蕩的動(dòng)力學(xué)特性感興趣.Barron 和Sen[35]考慮了在一質(zhì)量-彈簧結(jié)構(gòu)上4 個(gè)耦合van der Pol 振蕩的同步性.Siewe 等[36]從實(shí)驗(yàn)上分析了在非線性與時(shí)滯耦合的兩個(gè)van der Pol 振蕩的同步性,實(shí)驗(yàn)得到的結(jié)果與理論分析相一致.現(xiàn)考察兩自由度耦合van der Pol 振蕩方程[37-38]

其 中,μ1=-0.1 ,μ2=0.5 ,γ1=0.1 ,γ2=0.08 .同樣 利用IHB 法求解這類強(qiáng)非線性耦合振蕩系統(tǒng),引入新的時(shí)間變量 τ =ωt,方程變?yōu)?/p>

因方程(39)只含有奇次項(xiàng),系統(tǒng)是對(duì)稱的振蕩響應(yīng),可將x1和x2展開為含有m個(gè)奇次諧波項(xiàng)的傅里葉級(jí)數(shù)

其中,m=10 .根據(jù)式(14),同樣可得方程(36),其中,KA是 4m×4m的矩陣,對(duì)應(yīng)于傅里葉系數(shù)向量

同樣選定 Δ λ 為主動(dòng)增量,再選定式(41)中的諧波項(xiàng)系數(shù)b11為參考量,那么該參考量的增量為零,即 Δb11=0 ,根據(jù)方程(36)即可得在不同參數(shù) λ 下系統(tǒng)響應(yīng)的頻率 ω ,即 λ -ω 響應(yīng)曲線.在此例中,給定不同的初值,利用IHB 法會(huì)得到不同的 λ -ω 響應(yīng)曲線.設(shè)選定初值為ω0=1.0

其中,l為取1 的個(gè)數(shù),m-l表示取0 的個(gè)數(shù).圖6 表示耦合van der Pol 振蕩方程(38)在其系數(shù) μ1=-0.1,μ2=0.5,γ1=0.1 和 γ2=0.08 ,初值條件為 ω0=1.0 和式(42)下,3 種IHB 法計(jì)算所得的系統(tǒng)頻率 ω 與參數(shù) λ 的對(duì)應(yīng)關(guān)系曲線(ω -λ),在這些曲線上的每一點(diǎn)就對(duì)應(yīng)系統(tǒng)的一個(gè)極限環(huán).在圖6 中,3 種IHB 法都可求得系統(tǒng)響應(yīng)頻率 ω ≈1.0 附近的 ω -λ 關(guān)系曲線,兩種改進(jìn)的GIHB 法可求得系統(tǒng)響應(yīng)頻率ω ≈0.67 和ω ≈0.3 附近的的ω -λ 關(guān)系曲線,而結(jié)合狗腿法和BLS 算法的GIHB2 法還可求得系統(tǒng)響應(yīng)頻率 ω ≈0.4 和 ω ≈0.2 附近的的 ω -λ 關(guān)系曲線.表明了兩種改進(jìn)后的GIHB 法相比傳統(tǒng)的IHB 法有更廣的初值選定范圍,得到了傳統(tǒng)IHB 法較難求得的解,特別是結(jié)合狗腿法和BLS 算法的GIHB2 法可以得到系統(tǒng)解的全貌.值得一提的是:進(jìn)一步研究表明,在 ω ≈1.0 附近的 ω -λ 關(guān)系曲線有一半的點(diǎn)對(duì)應(yīng)的極限環(huán)是穩(wěn)定的,另一半對(duì)應(yīng)的極限環(huán)是不穩(wěn)定的;其余4 條曲線對(duì)應(yīng)的極限環(huán)也是不穩(wěn)定的.

圖6 不同初值條件下用3 種IHB 法求得耦合范德波爾振蕩方程(38)的解,其中,μ 1=-0.1 ,μ 2=0.5 ,γ 1=0.1 和γ2=0.08Fig.6 The solutions of coupled van der Pol Eq.(38) with different initial values using the three IHB methods,where μ 1=-0.1 ,μ 2=0.5,γ1=0.1 ,andγ2=0.08

圖7 所示為不同初值條件下用3 種IHB 法求解耦合van der Pol 振蕩方程(38)在 ω ≈1.0 附近迭代收斂情況,其中,λ =0.08 ,μ1=-0.1,μ2=0.5 ,γ1=0.1和 γ2=0.08 .圖7(a)中選定初值為頻率 ω0=1.0,諧波項(xiàng)系數(shù)取式(42)中的l=7 ,b11=0.7,可以看出三種方法都可收斂,傳統(tǒng)IHB 法誤差 ‖R‖ 在前10 迭代過程中在 [ 0,100] 間波動(dòng),然后開始下降,達(dá)到收斂條件.用兩種改進(jìn)的GIHB 法得到的誤差 ‖R‖ 一直減少,直至滿足迭代條件.看出傳統(tǒng)IHB 法需要18 步迭代達(dá)到收斂條件,而兩種改進(jìn)的GIHB 法分別只需要12 步和10 步迭代達(dá)到收斂條件.圖7(b)中選定初值為頻率 ω0=1.0 ,諧波項(xiàng)系數(shù)取式(42)中的l=8,b11=0.8,可以看出傳統(tǒng)IHB 法在迭代的前10 步過程中波動(dòng),第11 步誤差徒然變大,之后迭代不收斂.而兩種改進(jìn)的GIHB 法得到的誤差 ‖R‖ 一直減少,直至滿足迭代條件.改進(jìn)后的GIHB1 法和GIHB2 法分別只需要14 步和11 步迭代達(dá)到收斂條件.圖8所示為取λ =0.08 ,μ1=-0.1 ,μ2=0.5 ,γ1=0.1 和γ2=0.08時(shí)耦合van der Pol 振蕩方程(38)的極限環(huán),可以看出3 種IHB 法求得的結(jié)果一致,并與用四階R-K 法求得的數(shù)值結(jié)果高度吻合.

圖7 不同初值條件下用3 種IHB 法求解耦合van der Pol 振蕩方程(38)在 ω ≈1.0 附近迭代收斂情況,其中,λ =0.08 ,μ 1=-0.1,μ2=0.5,γ 1=0.1 和γ2=0.08Fig.7 Iteration convergence of solution of coupled van der Pol Eq.(38)near ω ≈1.0 with different initial values using the three IHB methods,where λ =0.08 ,μ 1=-0.1 ,μ 2=0.5 ,γ 1=0.1 ,andγ2=0.08

圖8 用3 種IHB 法和四階龍格-庫(kù)塔數(shù)值法求得耦合van der Pol 振蕩方程(38)的解,其中,λ =0.08 ,μ 1=-0.1 ,μ 2=0.5 ,γ 1=0.1 和γ2=0.08Fig.8 Solutions of coupled van der Pol Eq.(38) by the three IHB methods and numerical integration using the fourth-order R-K method,where λ =0.08 ,μ 1=-0.1 ,μ 2=0.5 ,γ 1=0.1,and γ2=0.08

5 結(jié)論

針對(duì)傳統(tǒng)IHB 法遇到的初值選擇問題,本文提出了兩種改進(jìn)的IHB 法,第一種是引入BLS 優(yōu)化算法的改進(jìn)IHB 法(GIHB1),其目的是為了調(diào)節(jié)IHB 法的迭代步長(zhǎng);第二種是在第一種改進(jìn)的基礎(chǔ)上再結(jié)合狗腿法思想的改進(jìn)IHB 法(GIHB2),其目的是為了調(diào)節(jié)迭代方式,使迭代方向沿著較快的下降方向,從而加快BLS 算法的迭代運(yùn)算速度.兩個(gè)自治系統(tǒng)的算例說明了兩種改進(jìn)IHB 法的有效性,第一算例表明兩種改進(jìn)的IHB 法的初值選擇范圍遠(yuǎn)大于傳統(tǒng)IHB 法的初值選擇范圍,且GIHB2 法在一些情況下求得解的收斂步數(shù)遠(yuǎn)少于GIHB1 法;第二個(gè)算例表明在給定的初值條件下,傳統(tǒng)的IHB 法只能得到1 條 ω -λ 關(guān)系曲線,GIHB1 法可以得到3 條ω-λ關(guān)系曲線,而結(jié)合狗腿法思想的GIHB2 法可以得到5 條 ω -λ 關(guān)系曲線,同樣地,GIHB2 法求得解的收斂步數(shù)也小于GIHB1 法.另外,兩種改進(jìn)的GIHB 法只是改進(jìn)了增量過程,而沒有改變諧波平衡過程,計(jì)算的結(jié)果與傳統(tǒng)IHB 的結(jié)果是一致的,所以兩種改進(jìn)的GIHB 法對(duì)于非自治系統(tǒng)的分析也是有效的.綜上,引入BLS 優(yōu)化算法再結(jié)合狗腿法思想的GIHB2 法在解決初值問題上是傳統(tǒng)IHB 法的有效補(bǔ)充,可進(jìn)一步推廣至高維非線性振動(dòng)系統(tǒng)上.