張利娟 張華彪 李欣業(yè)

1)(河北工業(yè)大學(xué)機(jī)械工程學(xué)院,天津 300130)

2)(天津商業(yè)大學(xué)機(jī)械工程學(xué)院,天津 300134)

(2018年9月8日收到;2018年10月22日收到修改稿)

針對(duì)基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺的非線性動(dòng)力學(xué)響應(yīng)進(jìn)行研究,利用拉格朗日方程建立了系統(tǒng)的動(dòng)力學(xué)方程.將離散傅里葉變換、諧波平衡法以及同倫延拓方法相結(jié)合,對(duì)系統(tǒng)的周期響應(yīng)進(jìn)行求解,避免了傳統(tǒng)方法計(jì)算中使用泰勒展開引起的小振幅的限制,與數(shù)值計(jì)算結(jié)果的對(duì)比表明該求解方法具有較高的精確度.利用Floquet理論分析了周期響應(yīng)的穩(wěn)定性,給出了基礎(chǔ)運(yùn)動(dòng)振幅和頻率對(duì)系統(tǒng)周期響應(yīng)的影響.研究發(fā)現(xiàn):對(duì)應(yīng)某些基礎(chǔ)頻率和振幅,系統(tǒng)的周期響應(yīng)可能發(fā)生Hopf分岔;利用數(shù)值計(jì)算研究了Hopf分岔后系統(tǒng)響應(yīng)隨基礎(chǔ)頻率和振幅的變化,發(fā)現(xiàn)系統(tǒng)出現(xiàn)了倍周期運(yùn)動(dòng)、擬周期運(yùn)動(dòng)和混沌等復(fù)雜的動(dòng)力學(xué)行為.研究表明系統(tǒng)進(jìn)入混沌的主要路徑是擬周期環(huán)面破裂和陣發(fā)性.

1 引 言

彈簧擺作為一個(gè)典型的兩自由度非線性系統(tǒng),在許多關(guān)于非線性振動(dòng)的專著中都有所論述[1,2],它既可以看成是單擺的推廣,也可以看成是單自由度彈簧質(zhì)量系統(tǒng)的推廣.由于人們對(duì)這兩類單自由度系統(tǒng)非常熟悉,所以彈簧擺特別適合于說明非線性振動(dòng)系統(tǒng)中的內(nèi)共振以及其他動(dòng)力學(xué)現(xiàn)象.事實(shí)上,在研究大型起重船或吊車作業(yè)時(shí)懸吊物的擺振問題時(shí)就可以將此簡化為彈簧擺模型,因此其非線性動(dòng)力學(xué)問題一直得到研究者的關(guān)注.

Eissa等[3]用多尺度法研究了六次非線性模型彈簧擺在兩個(gè)模態(tài)都受到簡諧激勵(lì)時(shí)主共振、亞諧共振、超諧共振、組合共振響應(yīng)的四階近似解,發(fā)現(xiàn)參數(shù)的變化會(huì)導(dǎo)致頻響曲線的側(cè)彎.Alasty和Shabani[4]研究了1:2內(nèi)共振彈簧擺沿呼吸和擺動(dòng)模態(tài)都受簡諧激勵(lì)時(shí)主共振響應(yīng)的混沌運(yùn)動(dòng),發(fā)現(xiàn)除了奇怪吸引子外還有多個(gè)共存的規(guī)則吸引子,且后者的吸引域邊界是分形的.Starosta等[5]研究了內(nèi)共振彈簧擺受不同頻率簡諧激勵(lì)時(shí)的參數(shù)激勵(lì)共振-主共振響應(yīng),結(jié)果表明頻響曲線既可以呈硬特性也可能呈軟特性.Awrejcewicz等[6]對(duì)擺錘沿?cái)[動(dòng)的徑向及切向分別受簡諧激勵(lì)的彈簧擺,利用多尺度法研究了彈簧的非線性和耦合強(qiáng)度對(duì)穩(wěn)態(tài)和瞬態(tài)共振響應(yīng)的影響,得到了穩(wěn)態(tài)共振響應(yīng)的頻響特性方程以及非穩(wěn)態(tài)響應(yīng)的近似漸近解.Klimenko等[7]研究了彈簧擺和一個(gè)單擺吸振器系統(tǒng)的動(dòng)力學(xué)響應(yīng),計(jì)算了系統(tǒng)的非線性模態(tài),并通過穩(wěn)定性分析給出了不同非線性模態(tài)的穩(wěn)定和不穩(wěn)定區(qū)域.蕭寒等[8]基于多尺度法研究了在單頻激勵(lì)下1:2內(nèi)共振平方非線性模型彈簧擺主共振響應(yīng)鞍結(jié)分岔的控制問題,設(shè)計(jì)的反饋控制器不僅能使系統(tǒng)的振幅得到有效的控制,還可以將彈簧擺的擾動(dòng)控制到一條穩(wěn)定的周期運(yùn)動(dòng)軌道上.全局響應(yīng)方面,Sousa等[9]提出了一種用于研究非線性耦合系統(tǒng)能量交換方式以及系統(tǒng)參數(shù)對(duì)能量交換影響的方法,并將其應(yīng)用于彈簧擺系統(tǒng),分析了系統(tǒng)的擺動(dòng)和呼吸運(yùn)動(dòng)之間的能量轉(zhuǎn)換,給出了彈簧擺系統(tǒng)能量分布的全局特征.Lee[10]利用插值映射法研究了平方非線性模型彈簧擺存在1:2內(nèi)共振關(guān)系時(shí)在簡諧激勵(lì)下主共振穩(wěn)定周期響應(yīng)的全局吸引域問題,發(fā)現(xiàn)在四維狀態(tài)空間中存在一個(gè)特殊的平面,包含了全部吸引子.Lee和Park[11,12]利用多尺度法分別研究了同時(shí)含有平方和立方非線性并具有1:2內(nèi)共振關(guān)系的彈簧擺在簡諧激勵(lì)下主共振響應(yīng)的一階和二階近似解的混沌運(yùn)動(dòng),分岔圖和李雅普諾夫指數(shù)結(jié)果表明,二階近似解比一階近似解能更好地與原系統(tǒng)吻合.Zaki等[13]利用數(shù)值方法討論了只沿切線方向受簡諧激勵(lì)的彈簧擺非主共振響應(yīng)的分岔和通向混沌的道路,發(fā)現(xiàn)彈簧的非線性對(duì)其有重要的影響.Tian等[14,15]研究了具有無理式和分式形式非線性恢復(fù)力的彈簧擺,利用多尺度法分析了泰勒展開后的系統(tǒng)的周期響應(yīng),利用改進(jìn)的Melnikov方法給出了系統(tǒng)出現(xiàn)混沌的閾值,并進(jìn)行了數(shù)值驗(yàn)證.Awrejcewicz等[16]研究了同時(shí)有兩個(gè)方向外激勵(lì)的彈簧擺的非穩(wěn)態(tài)響應(yīng),利用多尺度法給出了系統(tǒng)非穩(wěn)態(tài)響應(yīng)的包絡(luò)線和極限相軌跡.Digilov等[17]建立了彈簧擺上變質(zhì)量阻尼振蕩模型,該模型的質(zhì)量以恒定速率衰減,討論了質(zhì)量損失率對(duì)阻尼參數(shù)大小的影響.Eissa等[18]將橫向調(diào)諧吸振器應(yīng)用于多參數(shù)激振的彈簧擺系統(tǒng)的振動(dòng)控制中,研究了次諧共振和內(nèi)共振情況下的周期解及其穩(wěn)定性,討論了吸振器和系統(tǒng)參數(shù)對(duì)響應(yīng)的影響.

除了簡諧激勵(lì)外,受基礎(chǔ)激勵(lì)的彈簧擺也是研究的熱點(diǎn)之一.Gitterman[19]對(duì)懸掛點(diǎn)做豎直方向簡諧運(yùn)動(dòng)和受到外加激勵(lì)的彈簧擺進(jìn)行了比較研究,為此將運(yùn)動(dòng)微分方程中的三角函數(shù)做泰勒展開后截取到四次項(xiàng),用多尺度法求取了四階近似解.結(jié)果表明,當(dāng)外激勵(lì)和黏性阻尼為同階小時(shí),二者的低階近似解是一樣的,但當(dāng)阻尼力比外激勵(lì)更小時(shí),低階近似解是不同的.Amer等[20,21]分別對(duì)沿圓形軌跡和橢圓軌跡移動(dòng)的彈簧擺進(jìn)行研究,采用多尺度法求得了系統(tǒng)的周期解,并進(jìn)行了穩(wěn)定性分析,同時(shí)利用數(shù)值計(jì)算分析了周期解失穩(wěn)后的響應(yīng).他們又將擺球由質(zhì)點(diǎn)改為剛體開展研究[22],對(duì)系統(tǒng)的共振情況進(jìn)行分類,給出了穩(wěn)態(tài)解的可解條件.

工程中很多彈簧擺系統(tǒng)在工作時(shí)其基礎(chǔ)是水平運(yùn)動(dòng)的,可能是在平衡位置附近的左右往復(fù)運(yùn)動(dòng),也可能是行進(jìn)過程中的變速運(yùn)動(dòng),然而針對(duì)基礎(chǔ)水平運(yùn)動(dòng)彈簧擺的研究并不多見.同時(shí)對(duì)于已有的彈簧擺甚至單擺的研究,特別是周期解的計(jì)算中,由于動(dòng)力學(xué)方程中存在擺角的三角函數(shù)項(xiàng),傳統(tǒng)解析方法無法處理,只能對(duì)擺角進(jìn)行泰勒展開后再進(jìn)行研究.眾所周知,泰勒展開只在展開點(diǎn)附近是相對(duì)精確的,因此泰勒展開后的動(dòng)力學(xué)方程只適用于擺動(dòng)角度很小的情況,而在擺動(dòng)角度很大時(shí),其所得的分析結(jié)果往往是不準(zhǔn)確,甚至錯(cuò)誤的.特別是對(duì)于基礎(chǔ)激勵(lì)的情況,其動(dòng)力學(xué)方程中除回復(fù)力外,在激勵(lì)項(xiàng)中也含有擺角的三角函數(shù),泰勒展開不僅改變了回復(fù)力的大小,還改變了激勵(lì)的形式.因此針對(duì)未經(jīng)泰勒展開的彈簧擺系統(tǒng)的動(dòng)力學(xué)方程開展研究,特別是研究大擺角情況下的周期解及其穩(wěn)定性具有更高的理論價(jià)值.

基于上述問題,本文針對(duì)基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺的響應(yīng)開展研究.第2部分利用拉格朗日方程建立了基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺的動(dòng)力學(xué)方程;第3部分將諧波平衡法、離散傅里葉變換以及同倫延拓方法相結(jié)合,直接對(duì)未經(jīng)泰勒展開彈簧擺系統(tǒng)的動(dòng)力學(xué)方程進(jìn)行求解,求得了系統(tǒng)的周期響應(yīng);第4部分針對(duì)系統(tǒng)周期解的穩(wěn)定性進(jìn)行研究;第5部分分析了基礎(chǔ)頻率和振幅對(duì)系統(tǒng)響應(yīng)的影響;最后對(duì)本文的主要結(jié)論進(jìn)行總結(jié).

2 基礎(chǔ)水平運(yùn)動(dòng)彈簧擺的動(dòng)力學(xué)模型

圖1 基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺Fig.1.Spring pendulum with base motion in the horizontal direction.

圖1 為基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺模型,其中x1為擺球沿彈簧長度方向的位移,x2為彈簧擺的角位移,m和k分別表示擺球的質(zhì)量和彈簧的剛度,l為擺長,有l(wèi)=l0+mg/k,其中l(wèi)0為彈簧的原長.我們定義o為坐標(biāo)系原點(diǎn),建立如圖1所示x-o-y坐標(biāo)系,則擺球的位置可表示為

因此擺球的速度可寫作:

圖2 彈簧擺響應(yīng)的相圖和頻譜 (a)Ab=0.3 m,ω=4 rad/s,初值(1.09,1.43,1.26,0.11),仿真時(shí)間τ=0—800π;(b)Ab=0.3 m,ω=2.22 rad/s,初值(0.25,0.93,?0.33,0.98),仿真時(shí)間τ=0—800πFig.2.Phase diagram and spectrum of spring pendulum:(a)Ab=0.3 m,ω=4 rad/s,initial value(1.09,1.43,1.26,0.11),simulation time τ =0–800π;(b)Ab=0.3 m,ω =2.22 rad/s,initial value(0.25,0.93,?0.33,0.98),simulation time τ =0–800π.

擺球的動(dòng)能、勢(shì)能、瑞利耗散函數(shù)分別為利用拉格朗日方程可得系統(tǒng)的運(yùn)動(dòng)方程,有

本文主要考慮基礎(chǔ)做簡諧運(yùn)動(dòng)的情況,因此有x3=Abcosωt.定義X1=x1/l,X2=x2,τ=ωt,對(duì)方程進(jìn)行無量綱化并化簡,有

其中

對(duì)方程(5)進(jìn)行數(shù)值求解,如無特殊說明,本文設(shè)定計(jì)算參數(shù)為m=0.5 kg,k=20 N/m,l=1 m,g=10 m/s2,c1=0.1 Ns/m,c2=0.1 Ns.

圖2給出了彈簧擺運(yùn)動(dòng)的相圖和頻譜,可以看到系統(tǒng)周期響應(yīng)的頻率成分非常復(fù)雜,采用傳統(tǒng)方法求1次或者幾次諧波解無法準(zhǔn)確地描述周期解的響應(yīng)特性,因此本文將諧波平衡法和同倫延拓方法相結(jié)合對(duì)系統(tǒng)進(jìn)行高次諧波求解.

3 基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺的周期解

設(shè)方程的解為

其中A和B分別為各次諧波的系數(shù),n為所求解的最高次諧波數(shù).將(7)式代入方程(5),利用傅里葉變換提取各次諧波系數(shù),可得關(guān)于振幅的代數(shù)方程如下:

下面將利用同倫延拓方法對(duì)(8)式進(jìn)行求解,為了便于說明求解過程,令z=[z1,z2,···,zh],其中h=4n+3,z1=A0,z2=B0,z4i?1=A1i,z4i=A2i,z4i+2=B2i,z4i+1=B1i,zh為所要研究的系統(tǒng)參數(shù),則方程(8)可寫作

其中,H(z)=[H1,H2,···,Hh?1]T. 設(shè)定zj(zj∈Rh)為解曲線上的一個(gè)正則點(diǎn)(z0可以通過給定系統(tǒng)參數(shù),利用牛頓拉普森法求得),沿該點(diǎn)的切線方向預(yù)估,下一個(gè)預(yù)測(cè)點(diǎn)可寫作

式中δj為第j部的步長;uj為預(yù)測(cè)方向,滿足下列要求

從p=1開始迭代,直到滿足∥yp+1?yp∥6ε時(shí)停止,下一個(gè)點(diǎn)zj+1=yp+1.

理論上利用上述方法可以求得系統(tǒng)參數(shù)對(duì)各次諧波振幅的影響.但實(shí)際計(jì)算中,H(yp)和H′(y0)的計(jì)算需要對(duì)(8)式中的積分進(jìn)行直接解析求解,考慮到f1,f2中含有X2的三角函數(shù)項(xiàng),這是非常困難的,特別是當(dāng)諧波數(shù)較多,即n較大時(shí),解析求解幾乎是不可能的.因此我們采用數(shù)值方法求取H(yp)和H′(y0),將τ在一個(gè)周期2π內(nèi)等分成q份,則τ變?yōu)橐粋€(gè)等差序列[0,2π/q,4π/q,6π/q,···,2(q?1)π/q],給定yp的值,代入到(7)式中,則X1,X2都變成了時(shí)間序列,然后將X1,X2代入方程(5),則f1,f2變?yōu)?h?1)×q的序列,對(duì)該序列進(jìn)行離散傅里葉變換可以求得H(yp)的值.對(duì)于導(dǎo)數(shù)矩陣H′(y0),有

矩陣中各項(xiàng)可以通過差分形式求得,有

利用上述方法對(duì)方程進(jìn)行求解時(shí),理論上n取值越大,所得解的精度越高,然而n取值過大,則會(huì)影響計(jì)算效率.我們對(duì)系統(tǒng)的周期解進(jìn)行了大量的數(shù)值計(jì)算,發(fā)現(xiàn)系統(tǒng)的周期解的主要頻率成分在20次以內(nèi),因此本文設(shè)解的諧波次數(shù)n=20.圖3給出了基礎(chǔ)激勵(lì)的彈簧擺系統(tǒng)的計(jì)算結(jié)果,其中實(shí)線為20次諧波解,圓點(diǎn)為龍格-庫塔方法的計(jì)算結(jié)果,可以看到兩者吻合得非常好,特別是圖3(b)中擺動(dòng)方向的振幅已經(jīng)達(dá)到了1.5(角度值約85.9?),遠(yuǎn)遠(yuǎn)超出了微幅振動(dòng)的范圍,因此本文的方法不僅具有較高的求解精度,而且適合于彈簧擺大幅振動(dòng)的情況.

圖3 20次諧波解和數(shù)值解的對(duì)比(HB,諧波平衡法的結(jié)果;RK,龍格-庫塔法結(jié)果) (a)ω=4 rad/s,Ab=0.3 m;(b)ω=2.22 rad/s,Ab=0.3 mFig.3.Comparison of 20th harmonic solution and numerical solution(HB,the solution of harmonic balance method;RK,the solution given by Runge-Kutta method):(a)ω=4 rad/s,Ab=0.3 m;(b)ω=2.22 rad/s,Ab=0.3 m.

4 周期解的穩(wěn)定性分析

下面利用Floquet理論研究周期解的穩(wěn)定性,將方程(5)寫作狀態(tài)方程形式:

其中M=DF().顯然M為時(shí)變參數(shù)矩陣,利用Hsu方法求解,以單位陣為初值,在一個(gè)周期內(nèi)進(jìn)行數(shù)值積分,可求得Floquet乘數(shù)矩陣,其特征值為Floquet乘子.當(dāng)Floquet乘子從1處穿出單位圓時(shí),系統(tǒng)發(fā)生靜態(tài)分岔,Floquet乘子從?1處穿出單位圓時(shí),系統(tǒng)發(fā)生倍周期分岔,Floquet乘子從其他位置穿出單位圓時(shí),系統(tǒng)發(fā)生Hopf分岔.

5 基礎(chǔ)運(yùn)動(dòng)頻率和振幅的影響分析

圖4 基礎(chǔ)運(yùn)動(dòng)頻率對(duì)系統(tǒng)振幅的影響(SP,靜態(tài)分岔點(diǎn);HP,Hopf分岔點(diǎn))(a)Ab=0.02 m;(b)Ab=0.05 m;(c)Ab=0.25 m;(d)Ab=0.3 m;(e)Ab=0.4 m;(f)Ab=0.6 mFig.4.Influence of base frequency on the amplitude of system(SP,the static bifurcation point;HP,the Hopf bifurcation point):(a)Ab=0.02 m;(b)Ab=0.05 m;(c)Ab=0.25 m;(d)Ab=0.3 m;(e)Ab=0.4 m;(f)Ab=0.6 m.

圖4 給出了對(duì)應(yīng)不同基礎(chǔ)幅值時(shí)基礎(chǔ)運(yùn)動(dòng)頻率對(duì)系統(tǒng)振幅的影響曲線,這里的振幅是20次諧波疊加后振動(dòng)位移的最大值,表1是對(duì)應(yīng)圖4中各分岔點(diǎn)的分岔形式和Floquet乘子的變化,可以看到隨著基礎(chǔ)頻率的變化,系統(tǒng)的振幅出現(xiàn)了兩個(gè)峰值,當(dāng)基礎(chǔ)幅值A(chǔ)b=0.02 m時(shí),系統(tǒng)的響應(yīng)在整個(gè)頻率范圍內(nèi)都是穩(wěn)定的周期解.隨著基礎(chǔ)幅值的增大,響應(yīng)的兩個(gè)峰值分別向兩個(gè)方向發(fā)生偏斜,兩者之間的距離越來越遠(yuǎn),同時(shí)系統(tǒng)出現(xiàn)了不穩(wěn)定周期解.圖4(b)—(f)的變化趨勢(shì)是一致的,隨頻率的變化先發(fā)生兩次靜態(tài)分岔,基礎(chǔ)頻率在兩靜態(tài)分岔點(diǎn)之間時(shí)系統(tǒng)具有兩個(gè)穩(wěn)定的周期解,即系統(tǒng)具有雙穩(wěn)態(tài)現(xiàn)象.然后隨著基礎(chǔ)頻率的增大,系統(tǒng)將發(fā)生兩次Hopf分岔.而當(dāng)ω>3.5 rad/s時(shí),圖4(b)的下一個(gè)靜態(tài)分岔點(diǎn)出現(xiàn)在轉(zhuǎn)折點(diǎn)處,而圖4(c)—(f)中下一個(gè)靜態(tài)分岔點(diǎn)出現(xiàn)在曲線的中段.圖4(c)和圖4(d)在ω=5 rad/s附近有穩(wěn)定的周期解,然后在ω=6 rad/s附近通過靜態(tài)分岔失穩(wěn),而圖4(e)中當(dāng)ω>4.8 rad/s,系統(tǒng)不存在穩(wěn)定的周期解,圖4(f)中系統(tǒng)在ω=5 rad/s附近沒有穩(wěn)定的周期解,然而在ω=6 rad/s附近又出現(xiàn)了穩(wěn)定的周期解.利用數(shù)值計(jì)算對(duì)圖4(c)—(f)中通過靜態(tài)分岔點(diǎn)P1后的動(dòng)力學(xué)行為進(jìn)行研究,發(fā)現(xiàn)圖4(c)中通過P1后系統(tǒng)的響應(yīng)發(fā)生跳躍,跳到下方的周期解軌道上,而圖4(d)—(f)通過P1后系統(tǒng)的響應(yīng)不再是周期運(yùn)動(dòng),穩(wěn)定的運(yùn)動(dòng)中擺動(dòng)方向的位移隨著時(shí)間變化持續(xù)增大,如圖5所示.此時(shí)彈簧擺的運(yùn)動(dòng)表現(xiàn)為一種旋轉(zhuǎn)運(yùn)動(dòng),稱之為旋轉(zhuǎn)解,同時(shí)還可以看到此時(shí)呼吸運(yùn)動(dòng)的最大無量綱振幅超過了1010,由此可知旋轉(zhuǎn)解出現(xiàn)后彈簧擺必然會(huì)被破壞,這是工程中必須要避免的.數(shù)值計(jì)算發(fā)現(xiàn)在圖4(d)和圖4(f)中周期解通過P2點(diǎn)后也會(huì)變成旋轉(zhuǎn)解.

表1 對(duì)應(yīng)圖4的分岔類型和Floquet乘子的變化Table 1.Bifurcation type and change of Floquet multiplier corresponding to Fig.4.

圖5 旋轉(zhuǎn)解的時(shí)間歷程(Ab=0.3 m,ω=4.81 rad/s;初值(1.133,1.503,1.273,0.113);(c),(d),(e)為(a)的局部放大;(f),(g),(h)為(b)的局部放大)Fig.5.Time process of the rotational solution(Ab=0.3 m,ω=4.81 rad/s;initial value(1.133,1.503,1.273,0.113);panels(c),(d),(e)are the partial enlargement of panel(a),and panels(f),(g),(h)are the partial enlargement of panel(b)).

圖6 彈簧擺響應(yīng)隨基礎(chǔ)運(yùn)動(dòng)頻率變化的分岔圖和李雅普諾夫指數(shù)(Ab=0.3 m;ω=2.8 rad/s時(shí)的初值為(?0.27,0.34,?0.63,0.50);仿真時(shí)間τ=0—800π) (a)分岔圖;(b)李雅普諾夫指數(shù);(c),(d),(e),(f)分岔圖(a)的局部放大Fig.6.Bifurcation diagram and Lyapunov exponent of spring pendulum with the change of base frequency(Ab=0.3 m):(a)Bifurcation diagram;(b)Lyapunov exponent;(c),(d),(e),(f)enlargement of the bifurcation diagram(a).The initial value is(?0.27,0.34,?0.63,0.50)when ω =2.8 rad/s.Simulation time is τ=0–800π.

圖7 彈簧擺隨基礎(chǔ)運(yùn)動(dòng)頻率變化的典型動(dòng)力學(xué)響應(yīng)的龐加萊映射 (a)周期1運(yùn)動(dòng),ω=2.91 rad/s;(b)擬周期運(yùn)動(dòng),ω=2.92rad/s;(c)擬周期運(yùn)動(dòng),ω=2.93 rad/s;(d)混沌運(yùn)動(dòng),ω=2.9325 rad/s;(e)周期9運(yùn)動(dòng),ω=2.9346 rad/s;(f)混沌運(yùn)動(dòng),ω=2.9348 rad/s;(g)周期6運(yùn)動(dòng),ω=2.9632 rad/s;(h)混沌運(yùn)動(dòng),ω=2.9646 rad/s;(i)擬周期運(yùn)動(dòng),ω=2.992 rad/s;(j)混沌運(yùn)動(dòng),ω=3.002 rad/s;(k)周期7運(yùn)動(dòng),ω=3.17 rad/s;(l)混沌運(yùn)動(dòng),ω=3.172 rad/s;(m)周期27運(yùn)動(dòng),ω=3.1806 rad/s;(n)擬周期運(yùn)動(dòng),ω=3.182 rad/s;(o)擬周期運(yùn)動(dòng),ω=3.48 rad/sFig.7.Poincare map of typical dynamic behavior of spring pendulum with the change of base frequency:(a)Period 1 motion,ω=2.91 rad/s;(b)almost periodic motion,ω=2.92rad/s;(c)almost periodic motion,ω=2.93 rad/s;(d)chaos,ω=2.9325 rad/s;(e)period 9 motion,ω=2.9346 rad/s;(f)chaos,ω=2.9348 rad/s;(g)period 6 motion,ω=2.9632 rad/s;(h)chaos,ω=2.9646 rad/s;(i)almost periodic motion,ω=2.992 rad/s;(j)chaos,ω=3.002 rad/s;(k)period 7 motion,ω=3.17 rad/s;(l)chaos,ω=3.172 rad/s;(m)period 27 motion,ω=3.1806 rad/s;(n)almost periodic motion,ω =3.182 rad/s;(o)almost periodic motion,ω =3.48 rad/s.

對(duì)兩個(gè)Hopf分岔點(diǎn)之間的動(dòng)力學(xué)行為進(jìn)行數(shù)值計(jì)算,圖6給出了系統(tǒng)響應(yīng)隨基礎(chǔ)運(yùn)動(dòng)頻率變化的分岔圖和李雅普諾夫指數(shù),圖7給出了系統(tǒng)出現(xiàn)的典型動(dòng)力學(xué)響應(yīng)的龐加萊映射,可以看到系統(tǒng)在基礎(chǔ)運(yùn)動(dòng)頻率ω=2.91 rad/s時(shí),響應(yīng)為周期1運(yùn)動(dòng),隨著頻率的增大,周期1運(yùn)動(dòng)發(fā)生Hopf分岔變?yōu)閿M周期運(yùn)動(dòng)(見圖7(b)),頻率繼續(xù)增大后,擬周期運(yùn)動(dòng)發(fā)生環(huán)面倍化(見圖7c)),進(jìn)而擬周期環(huán)面破裂,系統(tǒng)的響應(yīng)進(jìn)入混沌(見圖7(d)).隨后系統(tǒng)的響應(yīng)在ω=2.933 rad/s附近退出混沌,變?yōu)橹芷?運(yùn)動(dòng)(見圖7(e)),然后在ω=2.935 rad/s附近通過陣發(fā)性再次進(jìn)入混沌.在ω=2.963 rad/s附近,系統(tǒng)退出混沌,響應(yīng)為周期6運(yùn)動(dòng),在ω=2.964 rad/s附近重新進(jìn)入混沌,此時(shí)系統(tǒng)的響應(yīng)在相平面上表現(xiàn)為6個(gè)小的吸引子,隨著頻率的增大,6個(gè)小的吸引子很快變成了一個(gè)大的吸引子.在ω=2.972 rad/s附近系統(tǒng)再次退出混沌,變?yōu)橹芷?運(yùn)動(dòng),在ω=2.984—2.994 rad/s時(shí)系統(tǒng)響應(yīng)為擬周期運(yùn)動(dòng),其龐加萊映射為9個(gè)不變環(huán),在ω=3.002 rad/s附近周期9運(yùn)動(dòng)通過陣發(fā)性回到混沌.在ω=3.168—3.215 rad/s之間,系統(tǒng)響應(yīng)在周期運(yùn)動(dòng)、混沌和擬周期運(yùn)動(dòng)之間發(fā)生多次切換,典型的響應(yīng)包括周期7運(yùn)動(dòng)、周期27運(yùn)動(dòng)、27不變環(huán)的擬周期運(yùn)動(dòng)和混沌(見圖7).然后系統(tǒng)在ω=3.4715 rad/s附近響應(yīng)變?yōu)榄h(huán)面倍化的擬周期運(yùn)動(dòng)(見圖7(o)),進(jìn)而回到一個(gè)不變環(huán)的擬周期運(yùn)動(dòng),隨著轉(zhuǎn)速的增大重新轉(zhuǎn)化為周期1運(yùn)動(dòng).需要注意的是,圖7(b)、圖7(c)、圖7(i)、圖7(n)和圖7(o)都是擬周期運(yùn)動(dòng),然而其運(yùn)動(dòng)的本質(zhì)是不同的,圖7(b)的擬周期運(yùn)動(dòng)中含有激勵(lì)的無量綱頻率1和Hopf分岔產(chǎn)生的不可約的頻率成分ωH及其線性組合,在圖7(c)中系統(tǒng)響應(yīng)發(fā)生了環(huán)面倍化,其實(shí)質(zhì)是ωH頻率成分發(fā)生了倍周期分岔,即響應(yīng)中出現(xiàn)了含ωH/2頻率的成分,圖7(o)的響應(yīng)形式和圖7(c)類似.而圖7(i)和圖7(n)的響應(yīng)實(shí)質(zhì)是由倍周期運(yùn)動(dòng)通過Hopf分岔產(chǎn)生的擬周期運(yùn)動(dòng),其頻率成分為1/N與Hopf分岔的頻率成分ωH及其線性組合,其中N是倍周期運(yùn)動(dòng)的周期數(shù).

我們同樣研究了基礎(chǔ)振幅對(duì)響應(yīng)的影響,圖8給出了對(duì)應(yīng)不同基礎(chǔ)運(yùn)動(dòng)頻率時(shí),基礎(chǔ)振幅對(duì)彈簧擺振幅的影響曲線,表2是與之對(duì)應(yīng)的各分岔點(diǎn)的分岔形式和Floquet乘子的變化.可以看到對(duì)應(yīng)于ω=2.3 rad/s和ω=4.6 rad/s兩種情況,曲線都出現(xiàn)了兩個(gè)靜態(tài)分岔點(diǎn),當(dāng)基礎(chǔ)振幅取某些值時(shí),系統(tǒng)出現(xiàn)了兩個(gè)穩(wěn)定解和一個(gè)不穩(wěn)定解,隨著基礎(chǔ)振幅的變化,響應(yīng)的振幅將發(fā)生跳躍現(xiàn)象;而在ω=6 rad/s時(shí),系統(tǒng)在Ab=0.2 m附近的靜態(tài)分岔點(diǎn)失穩(wěn)后將進(jìn)入旋轉(zhuǎn)狀態(tài),在Ab=0.53—0.6 m之間穩(wěn)定周期解和旋轉(zhuǎn)解共存,系統(tǒng)的響應(yīng)將取決于初始狀態(tài);在ω=2.9 rad/s時(shí),當(dāng)基礎(chǔ)振幅較小時(shí),系統(tǒng)響應(yīng)隨基礎(chǔ)振幅的變化將發(fā)生跳躍,而在Ab=0.328 m附近系統(tǒng)將發(fā)生Hopf分岔,響應(yīng)由周期運(yùn)動(dòng)變?yōu)閿M周期運(yùn)動(dòng).

圖8 基礎(chǔ)振幅對(duì)系統(tǒng)響應(yīng)振幅的影響(SP,靜態(tài)分岔點(diǎn);HP,Hopf分岔點(diǎn)) (a)ω=2.3 rad/s;(b)ω=2.9 rad/s;(c)ω=4.6 rad/s;(d)ω=6 rad/sFig.8.Influence of base amplitude on the amplitude of system(SP,the static bifurcation point;HP,the Hopf bifurcation point):(a)ω=2.3 rad/s;(b)ω=2.9 rad/s;(c)ω=4.6 rad/s;(d)ω=6 rad/s.

圖9 彈簧擺響應(yīng)隨基礎(chǔ)振幅變化的分岔圖和李雅普諾夫指數(shù)(ω=2.9 rad/s,Ab=0.3 m時(shí)的初值為(?0.30,0.34,?0.54,0.43),仿真時(shí)間τ=0—800π)Fig.9.Bifurcation diagram and Lyapunov exponent of spring pendulum with the change of base amplitude(ω=2.9 rad/s;the initial value is(?0.30,0.34,?0.54,0.43)when Ab=0.3 m;simulation time is τ=0–800π)

圖10 彈簧擺隨基礎(chǔ)振幅變化的典型動(dòng)力學(xué)行為的龐加萊映射 (a)周期1運(yùn)動(dòng),Ab=0.32 m;(b)擬周期運(yùn)動(dòng),Ab=0.34 m;(c)擬周期運(yùn)動(dòng),Ab=0.366 m;(d)混沌,Ab=0.368 m;(e)周期4運(yùn)動(dòng),Ab=0.372 m;(f)擬周期運(yùn)動(dòng),Ab=0.374 m;(g)周期16運(yùn)動(dòng),Ab=0.3752 m;(h)周期11運(yùn)動(dòng),Ab=0.38 m;(i)混沌,Ab=0.384 m;(j)擬周期運(yùn)動(dòng),Ab=0.4454 m;(k)周期7運(yùn)動(dòng),Ab=0.45 m;(l)混沌,Ab=0.485 mFig.10.Poincare map of typical dynamic behavior of spring pendulum with the change of base amplitude:(a)Period 1 motion,Ab=0.32 m;(b)almost periodic motion,Ab=0.34 m;(c)almost periodic motion,Ab=0.366 m;(d)chaos,Ab=0.368 m;(e)period 4 motion,Ab=0.372 m;(f)almost periodic motion,Ab=0.374 m;(g)period 16 motion,Ab=0.3752 m;(h)period 11 motion,Ab=0.38 m;(i)chaos,Ab=0.384 m;(j)almost periodic motion,Ab=0.4454 m;(k)period 7 motion,Ab=0.45 m;(l)chaos,Ab=0.485 m.

表2 圖8中的分岔類型和Floquet乘子的變化Table 2.Bifurcation type and change of Floquet multiplier corresponding to Fig.8.

為了了解系統(tǒng)響應(yīng)在Hopf分岔后隨基礎(chǔ)振幅的變化,我們利用龍格-庫塔方法對(duì)系統(tǒng)進(jìn)行數(shù)值求解.圖9給出了響應(yīng)隨基礎(chǔ)振幅變化的分岔圖和李雅普諾夫指數(shù)圖,可以看到系統(tǒng)響應(yīng)隨基礎(chǔ)振幅的變化出現(xiàn)了倍周期運(yùn)動(dòng)、擬周期運(yùn)動(dòng)和混沌等復(fù)雜的響應(yīng)形式.圖10給出了典型動(dòng)力學(xué)響應(yīng)的龐加萊映射,可以看到當(dāng)系統(tǒng)發(fā)生Hopf分岔后,響應(yīng)從周期1運(yùn)動(dòng)變?yōu)閿M周期運(yùn)動(dòng),隨著基礎(chǔ)振幅的增大,擬周期運(yùn)動(dòng)出現(xiàn)環(huán)面倍化(見圖10(c)),進(jìn)而環(huán)面發(fā)生破裂,系統(tǒng)響應(yīng)進(jìn)入混沌(見圖10(d)).隨著基礎(chǔ)振幅的進(jìn)一步增加,系統(tǒng)響應(yīng)退出混沌,變?yōu)橹芷?運(yùn)動(dòng)(見圖10(e)),然后通過Hopf分岔成為4個(gè)吸引不變環(huán)的擬周期運(yùn)動(dòng)(圖10(f))和周期16運(yùn)動(dòng)(圖10(h)).基礎(chǔ)振幅繼續(xù)增加,系統(tǒng)響應(yīng)出現(xiàn)了周期11運(yùn)動(dòng),然后通過陣發(fā)性進(jìn)入混沌(見圖10(i)).當(dāng)基礎(chǔ)振幅繼續(xù)增大時(shí),系統(tǒng)的響應(yīng)以混沌為主,只是在Ab=0.445—0.485 m之間,系統(tǒng)響應(yīng)在周期7運(yùn)動(dòng)、7個(gè)吸引不變環(huán)的擬周期運(yùn)動(dòng)和混沌之間出現(xiàn)了兩次切換.

綜合上面的分析,我們給出了基礎(chǔ)運(yùn)動(dòng)頻率-振幅平面上周期解的轉(zhuǎn)遷,可以看到靜態(tài)分岔邊界、Hopf分岔邊界和旋轉(zhuǎn)解邊界將參數(shù)平面分成了7個(gè)保持域,其中參數(shù)域②是擬周期運(yùn)動(dòng)、混沌以及復(fù)雜周期運(yùn)動(dòng)的分布域,參數(shù)域①內(nèi)系統(tǒng)是唯一的周期運(yùn)動(dòng),而參數(shù)域⑥內(nèi)為旋轉(zhuǎn)解.圖中陰影部分,即區(qū)域③,④,⑤,⑦,是初值敏感區(qū)域,即通常所說的雙穩(wěn)態(tài)區(qū)域,對(duì)應(yīng)不同的初值,系統(tǒng)可能收斂于不同的穩(wěn)態(tài)運(yùn)動(dòng).其中區(qū)域③和④對(duì)應(yīng)不同的初值,系統(tǒng)將收斂到振幅不同的周期解,而對(duì)于參數(shù)域⑤和⑦,當(dāng)初值較小時(shí),系統(tǒng)的響應(yīng)為振幅較小的周期運(yùn)動(dòng),當(dāng)初值較大時(shí),系統(tǒng)的響應(yīng)將會(huì)是連續(xù)的旋轉(zhuǎn)運(yùn)動(dòng).

圖11 基礎(chǔ)運(yùn)動(dòng)頻率-振幅平面上周期解的轉(zhuǎn)遷集(BS,靜態(tài)分岔邊界;HF,Hopf分岔邊界;BR,旋轉(zhuǎn)解邊界)Fig.11.Transition sets of the response on the base motion parameters plane(BS,boundary of static bifurcation;HF,boundary of Hopf bifurcation;BR,boundary of continuous rotational motion).

6 結(jié) 論

本文針對(duì)基礎(chǔ)水平運(yùn)動(dòng)的彈簧擺開展研究,利用拉格朗日方程建立了系統(tǒng)動(dòng)力學(xué)方程.將離散傅里葉變換、諧波平衡法以及同倫延拓方法相結(jié)合對(duì)系統(tǒng)的周期響應(yīng)進(jìn)行求解,突破了傳統(tǒng)方法使用泰勒展開導(dǎo)致的小振幅的限制,與數(shù)值計(jì)算結(jié)果的對(duì)比表明本文的求解方法具有較高的精確度.

研究了基礎(chǔ)運(yùn)動(dòng)的頻率和振幅對(duì)系統(tǒng)周期響應(yīng)的影響,發(fā)現(xiàn)基礎(chǔ)頻率對(duì)系統(tǒng)振幅的影響表現(xiàn)為雙峰,當(dāng)基礎(chǔ)振幅較大時(shí),在兩個(gè)峰值附近會(huì)發(fā)生跳躍現(xiàn)象.周期解振幅隨基礎(chǔ)振幅的增大而增大,對(duì)應(yīng)某些基礎(chǔ)頻率,周期解的振幅隨基礎(chǔ)振幅的變化也會(huì)發(fā)生跳躍.對(duì)應(yīng)某些基礎(chǔ)頻率和振幅,彈簧擺可能出現(xiàn)連續(xù)旋轉(zhuǎn)運(yùn)動(dòng),同時(shí)呼吸方向的振幅極大.基于上述結(jié)果,給出了基礎(chǔ)運(yùn)動(dòng)參數(shù)平面上響應(yīng)形式的轉(zhuǎn)遷.

對(duì)應(yīng)某些基礎(chǔ)頻率和振幅,系統(tǒng)的周期響應(yīng)可能發(fā)生Hopf分岔,利用數(shù)值仿真研究了Hopf分岔后響應(yīng)的變化,發(fā)現(xiàn)系統(tǒng)出現(xiàn)了倍周期運(yùn)動(dòng)、擬周期運(yùn)動(dòng)和混沌等復(fù)雜響應(yīng)形式.研究表明系統(tǒng)進(jìn)入混沌主要是通過擬周期環(huán)面破裂和陣發(fā)性.