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

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

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

?

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

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

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

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

強非線性,單擺,諧波-能量平衡法

引言

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

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

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

本文采用諧波-能量平衡法對強非線性單擺方程進行求解研究,并與KBM法進行了對比.

1　諧波-能量平衡法

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

研究形如

(1a)

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

(1b)

強非線性自由振動微分方程(1),如用一個等效的線性微分方程

(2a)

來代替.保持初始條件相同

(2b)

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

(2c)

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

(3)

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

(4)

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

(5a)

(5b)

(5c)

其中ψ=ωt.

1.1單項諧波-能量平衡法

設對稱性方程(1)的解為

x=a1cosψ+b1sinψ

(6)

用Ritz-Galerkin平均法：

(7a)

(7b)

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

(8)

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

1.2兩項諧波-能量平衡法

設對稱性方程方程(1)的解為：

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

(9)

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

(10a)

(10b)

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

(11)

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

2　單擺振動的周期解

2.1振動問題分類

單擺也稱為數(shù)學擺,其運動方程為

(12)

初始條件為：

(13)

其中固有角頻率和周期

(14)

(15)

ψ=ωt-φ

(16)

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

(17)

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

θ=acosψ

(18a)

其中：

(18b)

振動周期為

(18c)

(19)

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

ψ=ωt-φ

(20)

這是不需要考慮小參數(shù)的振動方程,稱為強非線性振動方程.

KBM法第一次近似解[1]

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

(21a)

幅—頻關系

(21b)

KBM法第二次近似解[1]

(22a)

幅—頻關系

(22b)

通常的振動問題按近似程度不同的工程要求可以分為：線性振動方程、弱非線性振動方程和強非線性振動方程.

圖1　單擺運動分類Fig. 1　Classification of the pendulum motion

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

2.2精確解

2.2.1單擺方程定性分析

(23)

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

(24)

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

(25)

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

圖2　單擺運動的相圖Fig. 2　The phase diagram of the pendulum motion

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

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

(26)

顯然,k<1對應于捕獲軌道,k>1對應于非捕獲軌道,k=1對應于界軌.

中心在(0,0),這是單擺振動的情況,設a是振幅,這時

(27)

而作用I的計算公式為

(28)

其中

(29)

引入變量α來代替q

(30)

表達式(28)可以改寫成

(31)

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

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

(32)

(33)

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

(34)

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

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

(35)

單擺振動的周期

(36)

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

(37)

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

(38)

又由(30)得

(39)

由(37)得

(40)

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

(41)

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

(42a)

(42b)

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

(43)

將式(36)與線性單擺運動的周期式(14)比較有

(44)

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

(45)

其中

(46)

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

(47)

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

(48)

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

(49)

2.3諧波-能量平衡法解

考察方程(16),設兩項諧波解為

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

(50)

令ψ=ωt用Ritz平均法：

(51)

① 單項諧波解幅—頻關系為

(52a)

初始條件的約束方程為：

(52b)

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

② 兩項諧波解幅—頻關系

(53a)

(53b)

初始條件的約束方程為：

(53c)

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

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

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

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

表1給出了諧波能量平衡法與其它方法的數(shù)值結(jié)果的比較.單項諧波法、兩項諧波法和精確解的振幅完全相同θmax=a,KBM漸近法對方程(22)的二次近似解振幅為θmax2,Karman對精確解的兩項級數(shù)展開法(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時單項諧波法、兩項諧波法與精確解的相圖比較.

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

單項諧波能量平衡法解：

θ=2.8cos(0.54097ω0t),

兩項諧波能量平衡法解：

θ=3.0236cos(0.52058ω0t)-

0.22361cos(1.5617ω0t),

KBM漸近法二階近似解：

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

Karman級數(shù)展開法二項近似解：

θ=3.0246cos(0.5610ω0t)-

0.25231cos(1.6830ω0t)

將近似解與精確解按角頻率比較,可知當a在160°附近時,單項諧波解與精確解的誤差為9.320%,兩項諧波解與精確解的誤差為5.200%,KBM法第二次近似解與精確解的誤差為15.17%,Karman法二項近似解與精確解的誤差為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)本文應用諧波-能量平衡法求解了強非線性單擺問題,分別給出了單項諧波解和兩項諧波解與精確解的比較.由圖4～圖9可見,當a>57°時,屬于強非線性,單項諧波解與精確解定性拓撲一致,而兩項諧波解與精確解相當一致.這表明諧波-能量平衡法既簡單,又精確.

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

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

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

1陳予恕.非線性振動.北京：高等教育出版社,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戴世強,莊峰青.一類非線性振動系統(tǒng)的漸近解.中國科學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李驪.強非線性振動系統(tǒng)的定性理論與定量方法.北京：科學出版社,1997 (Li L. Qualitative theory and quantitative method for strongly nonlinear vibration system. Beijing:Science Press, 1997 (in Chinese))

5陳樹輝,劉世齡,張佑啟,徐兆.強非線性振動的定量方法.廣州:廣東科技出版社,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李銀山,郝黎明,樹學鋒.強非線性Duffing方程的攝動解.太原理工大學學報,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è),劉波.分岔混沌非線性振動及其在工程中的應用.河北工業(yè)大學學報,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李銀山,陳予恕,吳志強.正交各向異性圓板非線性振動的亞諧分岔.機械強度, 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亞諧解.振動與沖擊,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李銀山,張善元,董青田,曹俊靈.用兩項諧波法求解強非線性Duffing方程.太原理工大學學報,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è),羅利軍.強非線性動力系統(tǒng)的兩項諧波法.太原理工大學學報, 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李銀山,張善元,劉波等.各種板邊條件下大撓度圓板自由振動的分岔解.機械強度, 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è).非對稱強非線性振動特征分析.動力學與控制學報,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收到修改稿.

*國家自然科學基金資助項目(10872063)