李銀山　李彤　韋炳威　李欣業(yè)

(1. 河北工業(yè)大學(xué)機(jī)械工程學(xué)院力學(xué)系, 天津　300130) (2. 華東理工大學(xué)承壓系統(tǒng)與安全教育部重點(diǎn)實(shí)驗(yàn)室, 上?！?00237)

?

用諧波-能量平衡法求解單擺方程?

李銀山1*李彤2韋炳威1李欣業(yè)1

(1. 河北工業(yè)大學(xué)機(jī)械工程學(xué)院力學(xué)系, 天津300130) (2. 華東理工大學(xué)承壓系統(tǒng)與安全教育部重點(diǎn)實(shí)驗(yàn)室, 上海200237)

應(yīng)用諧波-能量平衡法求解了強(qiáng)非線性單擺方程,諧波-能量平衡法與經(jīng)典的攝動(dòng)法和諧波平衡法不同,不是把微分方程和初始條件分離處理;而是把微分方程和初始條件同時(shí)處理.用諧波平衡,將描述動(dòng)力系統(tǒng)的二階常微分方程,化為以角頻率、振幅為變量的非線性代數(shù)方程組,考慮能量平衡,構(gòu)成角頻率、振幅為變量的封閉方程組求得解析解.諧波-能量平衡法將諧波平衡與能量平衡相結(jié)合,克服了二者的缺點(diǎn)吸取了二者的優(yōu)點(diǎn).實(shí)例表明,諧波-能量平衡法方法簡(jiǎn)單,取較少諧波就可以達(dá)到較高的精度.

強(qiáng)非線性,單擺,諧波-能量平衡法

引言

角頻率是描述周期振動(dòng)的最主要因素,采用通常的攝動(dòng)法[1-2]不能求解強(qiáng)非線性振動(dòng)問(wèn)題.近三十多年來(lái),強(qiáng)非線性振動(dòng)研究所取得的一系列成果,其突破點(diǎn)一般最終都可導(dǎo)出振動(dòng)頻率的瞬變性.比如時(shí)間變換法,橢圓函數(shù)法,頻閃法,推廣L-P法,等效線性化方法,改進(jìn)的多尺度法 ,FFT快速Galerkin法,增量諧波平衡法和攝動(dòng)增量法等[3-6].

張琪昌[7]等將待定固有頻率法與規(guī)范性方法相結(jié)合研究強(qiáng)非線性振動(dòng)問(wèn)題的求解.

李銀山2005年提出了求解強(qiáng)非線性振動(dòng)問(wèn)題的諧波-能量平衡法[13],其關(guān)鍵是采用諧波平衡加能量平衡構(gòu)成封閉的非線性代數(shù)方程組進(jìn)行求解.文獻(xiàn)[14-15]研究了采用諧波-能量平衡法求解對(duì)稱(chēng)強(qiáng)非線性動(dòng)力系統(tǒng)問(wèn)題.文獻(xiàn)[16]研究了采用諧波-能量平衡法求解非對(duì)稱(chēng)強(qiáng)非線性動(dòng)力系統(tǒng)問(wèn)題.

本文采用諧波-能量平衡法對(duì)強(qiáng)非線性單擺方程進(jìn)行求解研究,并與KBM法進(jìn)行了對(duì)比.

1　諧波-能量平衡法

諧波-能量平衡法的基本思想是把非線性微分方程組的解,用等效的線性微分方程組的解來(lái)解析逼近.首先采用諧波平衡,得到以振幅,角頻率為未知數(shù)的不完備非線性代數(shù)方程組(方程數(shù)小于未知數(shù));然后利用能量守恒原理,增加關(guān)于初始條件、振幅,角頻率之間協(xié)調(diào)的補(bǔ)充方程,從而構(gòu)成了關(guān)于振幅,角頻率為未知數(shù)的完備非線性代數(shù)方程組;對(duì)這個(gè)非線性代數(shù)方程組進(jìn)行求解,就可以得到近似解析解.

研究形如

(1a)

的振動(dòng)系統(tǒng).這里,f(x)是其變量的非線性奇函數(shù).初始條件為：

(1b)

強(qiáng)非線性自由振動(dòng)微分方程(1),如用一個(gè)等效的線性微分方程

(2a)

來(lái)代替.保持初始條件相同

(2b)

保持軌道的周期相同,能量相同.即

(2c)

H是系統(tǒng)的哈密頓能量函數(shù),h=const.設(shè)方程(1)的近似解析解為

(3)

將方程(1)式中的函數(shù)f(x)展開(kāi)成傅里葉級(jí)數(shù)

(4)

其中傅里葉系數(shù)為:

(5a)

(5b)

(5c)

其中ψ=ωt.

1.1單項(xiàng)諧波-能量平衡法

設(shè)對(duì)稱(chēng)性方程(1)的解為

x=a1cosψ+b1sinψ

(6)

用Ritz-Galerkin平均法：

(7a)

(7b)

根據(jù)能量平衡式(2c),得初始條件的約束方程為

(8)

由(7),(8)聯(lián)立可解得ω,a1,b1.

1.2兩項(xiàng)諧波-能量平衡法

設(shè)對(duì)稱(chēng)性方程方程(1)的解為：

x=a1cosψ+b1sinψ+a3cosψ+b3sinψ

(9)

令ψ=ωt用Ritz平均法便有：

(10a)

(10b)

根據(jù)能量平衡式(2c),得初始條件的約束方程為：

(11)

由(10),(11)聯(lián)立可解得ω,a1,b1,a3,b3.

2　單擺振動(dòng)的周期解

2.1振動(dòng)問(wèn)題分類(lèi)

單擺也稱(chēng)為數(shù)學(xué)擺,其運(yùn)動(dòng)方程為

(12)

初始條件為：

(13)

其中固有角頻率和周期

(14)

(15)

ψ=ωt-φ

(16)

其中Jk(a)為貝塞爾函數(shù).

(17)

這是線性振動(dòng)方程.方程(17)的解為

θ=acosψ

(18a)

其中：

(18b)

振動(dòng)周期為

(18c)

(19)

這是弱非線性振動(dòng)方程(一般要求0<ε?1,這里ε=1/6),即著名的軟彈簧Duffing方程.

ψ=ωt-φ

(20)

這是不需要考慮小參數(shù)的振動(dòng)方程,稱(chēng)為強(qiáng)非線性振動(dòng)方程.

KBM法第一次近似解[1]

θ1=acosψ,ψ=ω1t-φ

(21a)

幅—頻關(guān)系

(21b)

KBM法第二次近似解[1]

(22a)

幅—頻關(guān)系

(22b)

通常的振動(dòng)問(wèn)題按近似程度不同的工程要求可以分為：線性振動(dòng)方程、弱非線性振動(dòng)方程和強(qiáng)非線性振動(dòng)方程.

圖1　單擺運(yùn)動(dòng)分類(lèi)Fig. 1　Classification of the pendulum motion

由上分析可知：?jiǎn)螖[運(yùn)動(dòng)由線性變?yōu)榉蔷€性,其運(yùn)動(dòng)形態(tài)由單一的周期運(yùn)動(dòng)變?yōu)槎鄻踊倪\(yùn)動(dòng)了.

2.2精確解

2.2.1單擺方程定性分析

(23)

哈密頓函數(shù)為(取最低點(diǎn)為零勢(shì)能點(diǎn))

(24)

其中h=const為積分常數(shù).方程(23)就可寫(xiě)成

(25)

因此,非線性單擺系統(tǒng)(12)是一個(gè)保守系統(tǒng)或Hamilton系統(tǒng).qn=nπ是系統(tǒng)的平衡點(diǎn)(n為整數(shù)).當(dāng)n為偶數(shù)時(shí),qn為橢圓型不動(dòng)點(diǎn)坐標(biāo);當(dāng)n為奇數(shù)時(shí),qn為雙曲型不動(dòng)點(diǎn)坐標(biāo).從實(shí)際情況看,這樣的平衡位置只有兩個(gè)：一個(gè)是(0,0),中心點(diǎn),若給它以微小的位移,單擺作周期振蕩,平衡位置是穩(wěn)定的,它是單擺下垂,擺球位于下方的位置.另一個(gè)是(±π,0),鞍點(diǎn),若給它以微小的位移,單擺不再在平衡位置附近振蕩,而是旋轉(zhuǎn)起來(lái),平衡位置是不穩(wěn)定的,這是單擺擺球位于最上方的位置.

圖2　單擺運(yùn)動(dòng)的相圖Fig. 2　The phase diagram of the pendulum motion

2.2.2捕獲軌道、非捕獲軌道和界軌的解

為了給出單擺方程(12)一般相軌道的運(yùn)動(dòng)解,引進(jìn)能量參數(shù)

(26)

顯然,k<1對(duì)應(yīng)于捕獲軌道,k>1對(duì)應(yīng)于非捕獲軌道,k=1對(duì)應(yīng)于界軌.

中心在(0,0),這是單擺振動(dòng)的情況,設(shè)a是振幅,這時(shí)

(27)

而作用I的計(jì)算公式為

(28)

其中

(29)

引入變量α來(lái)代替q

(30)

表達(dá)式(28)可以改寫(xiě)成

(31)

其中K(k)和E(k)分別是第一類(lèi)和第二類(lèi)完全橢圓積分.

等式(31)確定了k的函數(shù)I.將它兩邊對(duì)k微分得

(32)

(33)

顯然,哈密頓函數(shù)H只依賴于I,由(26)和(31)確定,得到

(34)

其中k=k(I)是I=I(k)的反函數(shù),由(31)確定.

由(33)和(34)求得單擺振動(dòng)角頻率

(35)

單擺振動(dòng)的周期

(36)

正則變換q,p→ψ,I的母函數(shù),在變量替換(31)下為

(37)

其中F(α,k)和E(α,k)分別是第一類(lèi)和第二類(lèi)橢圓積分,α由等式(30)確定,而k=k(I)由(31)確定.角變量為

(38)

又由(30)得

(39)

由(37)得

(40)

考慮到(35)和(40),由公式(38)得出

(41)

由(29), (30)和(41)可得單擺振動(dòng)情況下引入作用—角變量的正則變換

(42a)

(42b)

其中sn()和cn()分別為Jacobi橢圓正弦和余弦函數(shù),且方程(12)的精確周期解為

(43)

將式(36)與線性單擺運(yùn)動(dòng)的周期式(14)比較有

(44)

將精確解(43)的右端[馮·卡門(mén)(Karman T V,1881～1963)]展開(kāi)成Fourier級(jí)數(shù),可得

(45)

其中

(46)

這是單擺旋轉(zhuǎn)的情況,方程(12)的精確解為

(47)

單擺旋轉(zhuǎn)的周期

(48)

這是單擺的同宿軌道.方程(12)的精確解為

(49)

2.3諧波-能量平衡法解

考察方程(16),設(shè)兩項(xiàng)諧波解為

θ=a1cos(ωt)+a3cos(3ωt)

(50)

令ψ=ωt用Ritz平均法：

(51)

① 單項(xiàng)諧波解幅—頻關(guān)系為

(52a)

初始條件的約束方程為：

(52b)

由(52)聯(lián)立可解得ω,a1.

② 兩項(xiàng)諧波解幅—頻關(guān)系

(53a)

(53b)

初始條件的約束方程為：

(53c)

由(53)聯(lián)立可解得ω,a1,a3.

2.4數(shù)值結(jié)果

圖3給出了諧波-能量平衡法與精確解周期隨振幅變化的關(guān)系T/T0～θmax的對(duì)比.

圖3　周期變化關(guān)系(°°°本法;——精確解)Fig. 3　Period-amplitude relationship(°°°:the present method's results; —— :the exact results)

表1給出了諧波能量平衡法與其它方法的數(shù)值結(jié)果的比較.單項(xiàng)諧波法、兩項(xiàng)諧波法和精確解的振幅完全相同θmax=a,KBM漸近法對(duì)方程(22)的二次近似解振幅為θmax2,Karman對(duì)精確解的兩項(xiàng)級(jí)數(shù)展開(kāi)法(45)振幅為θmaxii.

表1　諧波-能量平衡法與其它方法的數(shù)值結(jié)果比較

圖4　相圖a=π/3(°°°本法,——精確解)Fig. 4　Phase diagram a=π/3(°°°:the present method's results; —— :the exact results)

圖5　相圖a=π/2(°°°本法,——精確解)Fig. 5　phase diagram a=π/2(°°°:the present method's results; —— :the exact results)

圖6　相圖a=2π/3(°°°本法,——精確解)Fig. 6　Phase diagram a=2π/3(°°°:the present method's results; —— :the exact results)

圖4～圖9分別給出了a=π/3,a=π/2,a=2π/3,a=5π/6,a=2.8和a=3時(shí)單項(xiàng)諧波法、兩項(xiàng)諧波法與精確解的相圖比較.

以a=2.8為例,精確解為：

單項(xiàng)諧波能量平衡法解：

θ=2.8cos(0.54097ω0t),

兩項(xiàng)諧波能量平衡法解：

θ=3.0236cos(0.52058ω0t)-

0.22361cos(1.5617ω0t),

KBM漸近法二階近似解：

θ=2.8cos(0.5699ω0t)-0.242cos(1.7097ω0t),

Karman級(jí)數(shù)展開(kāi)法二項(xiàng)近似解：

θ=3.0246cos(0.5610ω0t)-

0.25231cos(1.6830ω0t)

將近似解與精確解按角頻率比較,可知當(dāng)a在160°附近時(shí),單項(xiàng)諧波解與精確解的誤差為9.320%,兩項(xiàng)諧波解與精確解的誤差為5.200%,KBM法第二次近似解與精確解的誤差為15.17%,Karman法二項(xiàng)近似解與精確解的誤差為13.37%.

圖7　相圖a=5π/6(°°°本法,——精確解)Fig. 7　Phase diagram a=5π/6(°°°:the present method's results; —— :the exact results)

圖8　相圖a=2.8(°°°本法,——精確解)Fig. 8　Phase diagram a=2.8(°°°:the present method's results; —— :the exact results)

圖9　相圖a=3(°°°本法,——精確解)Fig. 9　Phase diagram a=3(°°°:the present method's results; —— :the exact results)

3　結(jié)論

1)本文應(yīng)用諧波-能量平衡法求解了強(qiáng)非線性單擺問(wèn)題,分別給出了單項(xiàng)諧波解和兩項(xiàng)諧波解與精確解的比較.由圖4～圖9可見(jiàn),當(dāng)a>57°時(shí),屬于強(qiáng)非線性,單項(xiàng)諧波解與精確解定性拓?fù)湟恢?而兩項(xiàng)諧波解與精確解相當(dāng)一致.這表明諧波-能量平衡法既簡(jiǎn)單,又精確.

2)表1給出了諧波-能量平衡法、KBM漸近法、Karman級(jí)數(shù)展開(kāi)法和精確解的比較,結(jié)果表明,對(duì)a≤57°的弱非線性情況,KBM漸近法、Karman級(jí)數(shù)展開(kāi)法可以得到很好的結(jié)果.對(duì)a>57°的強(qiáng)非線性情況,KBM漸近法、Karman級(jí)數(shù)展開(kāi)法與精確解的偏離都比較大,而兩項(xiàng)諧波解與精確解相當(dāng)一致.這表明諧波-能量平衡法對(duì)強(qiáng)非線性系統(tǒng)可以得到很高的精度.

3)諧波-能量平衡法引入了能量平衡得出的初始條件約束方程,僅用兩項(xiàng)諧波就可得到較高的精度.克服了傳統(tǒng)的諧波平衡法需要取比較多的諧波數(shù)量才能得到較高精度的缺點(diǎn).

4)諧波-能量平衡法考慮了非線性等效特征.克服了傳統(tǒng)的等效線性化方法精度較差的缺點(diǎn).

1陳予恕.非線性振動(dòng).北京：高等教育出版社,2002 (Chen Y S. Nonlinear Vibration. Beijing: Higher education press, 2002 (in Chinese))

2Nayfeh A H, Mook D T. Nonlinear oscillations. New York: Wiley 1979

3戴世強(qiáng),莊峰青.一類(lèi)非線性振動(dòng)系統(tǒng)的漸近解.中國(guó)科學(xué)A輯,1986,29(1):34～40 (Dai S Q, Zhuang F Q. Asymptotic solution for a class of nonlinear oscillatory systems.TheScienceofChinaA, 1986,29(1):34～40 (in Chinese))

4李驪.強(qiáng)非線性振動(dòng)系統(tǒng)的定性理論與定量方法.北京：科學(xué)出版社,1997 (Li L. Qualitative theory and quantitative method for strongly nonlinear vibration system. Beijing:Science Press, 1997 (in Chinese))

5陳樹(shù)輝,劉世齡,張佑啟,徐兆.強(qiáng)非線性振動(dòng)的定量方法.廣州:廣東科技出版社,1997 (Chen S H, Liu S L, Zhang Y Q, Xu Z. Quantitative method for strong nonlinear vibration. Guangzhou: Guangdong Science and Technology Press, 1997 (in Chinese))

6李銀山,郝黎明,樹(shù)學(xué)鋒.強(qiáng)非線性Duffing方程的攝動(dòng)解.太原理工大學(xué)學(xué)報(bào),2000,31(5):516～520 (Li Y S, Hao L M, Shu X F. Asymptotic solution of strongly nonlinear Duffing equation.JournalofTaiyuanUniversityofTechnology, 2000,31(5):516～520 (in Chinese))

7Leung A Y T,Zhang Q C. Complex normal form for strongly nonlinear vibration systems exemplified by Duffing-van der Pol equation.JournalofSoundandVibration,1998,213(5):907～914

8李銀山,李欣業(yè),劉波.分岔混沌非線性振動(dòng)及其在工程中的應(yīng)用.河北工業(yè)大學(xué)學(xué)報(bào),2004,32(3):80～83 (Li Y S, Li X Y, Liu B. Bifurcation chaos nonlinear oscillations and their application in engineer.JournalofHebeiUniversityofTechnology, 2004,32(2):96～103 (in Chinese))

9李銀山,陳予恕,吳志強(qiáng).正交各向異性圓板非線性振動(dòng)的亞諧分岔.機(jī)械強(qiáng)度, 2001,23(2):148～151 (Li Y S, Chen Y S, Wu Z Q. Subharmonic bifurcation of nonlinear vibration of orthotropic circular plates.JournalofMechanicalStrength, 2001,23(2):148～151 (in Chinese))

11Li Y S, Zhang N M,Yang G T. 1/3 Subharmonic solution of elliptical sandwich plates.AppliedMathematicsandMechanics, 2003,24(10):1147～1157

12李銀山,張善元,張明路等.材料非線性圓板的1/2?1/4亞諧解.振動(dòng)與沖擊,2006,25(3):115～120 (Li Y S, Zhang S Y, Zhang M L. 1/2+1/4 Subharmonic solution of a circular plate from nonlinear of material.JournalofVibrationandShock, 2006,25(3):115～120 (in Chinese))

13李銀山,張善元,董青田,曹俊靈.用兩項(xiàng)諧波法求解強(qiáng)非線性Duffing方程.太原理工大學(xué)學(xué)報(bào),2005,36(6):690～693 (Li Y S, Zhang S Y, Dong Q T, Cao J L. Two harmonics method for strongly nonlinear Duffing equation.JounalofTaiyuanUniversityofTechnology, 2005,36(6):690～693 (in Chinese))

14李銀山,張善元,李欣業(yè),羅利軍.強(qiáng)非線性動(dòng)力系統(tǒng)的兩項(xiàng)諧波法.太原理工大學(xué)學(xué)報(bào), 2005,36(6):694～696 (Li YS, Zhang S Y, Li X Y, Luo L J. Two-harmonic method for strongly nonlinear dynamic systems.JournalofTaiyuanUniversityofTechnology, 2005,36(6):694～696 (in Chinese))

15李銀山,張善元,劉波等.各種板邊條件下大撓度圓板自由振動(dòng)的分岔解.機(jī)械強(qiáng)度, 2007,29(1):30～353 (Li Y S, Zhang S Y, Liu B. Bifurcate solutions of free vibration of a circular plate under various boundary conditions.JournalofMechanicalStrength, 2007,29(1):30～35 (in Chinese))

16李銀山,潘文波,吳艷艷,李欣業(yè).非對(duì)稱(chēng)強(qiáng)非線性振動(dòng)特征分析.動(dòng)力學(xué)與控制學(xué)報(bào),2012,10(1):15～20 (Li Y S, Pan W B, Wu Y Y, Li X Y. Asymmetric strongly nonlinear oscillation characteristic analysis.JournalofDynamicsandControl, 2012,10(1):15～20 (in Chinese))

*The project supported by the National Natural Science Foundation of China(10872063).

? Corresponding author E-mail: liyinshan@eyou.com

01 May 2015,revised 08 June 2015.

HARMONIC-ENERGY BALANCE METHOD FOR SOLVING PENDULUM EQUATION?

Li Yinshan1?Li Tong2Wei Bingwei1Li Xinye1

(1.DepartmentofMechanics,CollegeofMechanicalEng.,HebeiUniversityofTechnology,Tianjin300130China)(2.KeyLaboratoryofPressureSystemsandSafety,MinistryofEducation,EastChinaUniversityofScienceandTechnology,Shanghai200237,China)

A harmonic-energy balance method is put forward to solve the strong nonlinear pendulum equation. The difference from the classical perturbation method is that the harmonic balance method does not account the differential equation and initial conditions separately, but it considers both simultaneously. Through the harmonic-balance method, two-order ordinary differential equations describing dynamic systems become a set of nonlinear algebraic equations with the variables of angular frequency and amplitude. Considering the balance of energy, the close equations with angular frequency and amplitude as the variables can be solved. The harmonic-energy balance method is a combination of harmonic-balance and energy balance. It overcomes the shortcomings of both methods and takes their advantages. A case study also shows that the harmonic-energy balance method is simpler with higher precision although it takes less harmonics.

strong nonlinear,pendulum,harmonic-energy balance method

E-mail: liyinshan@eyou.com

10.6052/1672-6553-2015-047

2015-05-01收到第1稿,2015-06-08收到修改稿.

*國(guó)家自然科學(xué)基金資助項(xiàng)目(10872063)