Duffing振子倍周期分岔譜特性

唐元璋， 樓京俊， 翁雪濤, 朱石堅(jiān)

(1.海軍工程大學(xué) 動力工程學(xué)院，武漢 430033；2.船舶振動噪聲重點(diǎn)實(shí)驗(yàn)室，武漢 430033)

Feigenbaum[1-4]在所得非線性連續(xù)時(shí)間動力系統(tǒng)普適理論中對Duffing系統(tǒng)的譜結(jié)構(gòu)進(jìn)行過詳細(xì)研究，由第n次倍周期分岔開始，n次所有譜線將成為n+1次倍周期分岔后的偶數(shù)項(xiàng)元素，奇數(shù)項(xiàng)為新增譜線；n較大時(shí)新增譜線平均強(qiáng)度與上次倍周期分岔后新增譜線平均強(qiáng)度相比為常數(shù)μ=0.1525，即新增線譜平均強(qiáng)度較上次新增線譜平均強(qiáng)度下降8.18 dB。Wang等[5]亦獲得類似遞推結(jié)論。為驗(yàn)證推導(dǎo)的正確性，F(xiàn)eigenbaum對Rayleigh-Benard湍流系統(tǒng)的倍周期分岔行為分析、實(shí)驗(yàn)進(jìn)行常數(shù)驗(yàn)證。Linsay[6]在電路系統(tǒng)中觀察到電壓幅值倍周期分岔現(xiàn)象，已驗(yàn)證第3、4次倍周期分岔新增線譜平均強(qiáng)度分別接近常數(shù)μ。王本仁等[7]進(jìn)行液氮及水中次諧波實(shí)驗(yàn)，觀察到次諧波分岔現(xiàn)象。受噪聲影響及檢測設(shè)備限制，第5次分岔新增線譜較難觀察，故實(shí)驗(yàn)方法驗(yàn)證常數(shù)μ的精確性有限。倍周期分岔次數(shù)更多時(shí)，通常只用數(shù)值方法求解方程獲得功率譜，此需高頻分辨率譜分析才能觀察到微弱譜線。Tongue等[8-10]認(rèn)為數(shù)值方法有時(shí)因步長選取不當(dāng)會得到錯(cuò)誤結(jié)論。以上驗(yàn)證理論結(jié)論方法均為數(shù)值及實(shí)驗(yàn)方法，用近似解析方法求常數(shù)μ尚少見。增量諧波平衡法[11-13]為半數(shù)值半解析方法，可獲得非線性連續(xù)時(shí)間動力系統(tǒng)高精度高階諧波解[14]，既能克服實(shí)驗(yàn)方法精度難提高缺點(diǎn)，且無需數(shù)值方法求頻譜過程。本文用增量諧波平衡法以Duffing振子為研究對象，對常數(shù)μ進(jìn)行計(jì)算驗(yàn)證理論結(jié)論。其中理論推導(dǎo)以非線性連續(xù)時(shí)間動力系統(tǒng)展開，方法、過程可直接用于其它非線性連續(xù)時(shí)間動力系統(tǒng)。

1 倍周期分岔線譜變化規(guī)律

設(shè)N維非線性連續(xù)時(shí)間動力系統(tǒng)為：

(1)

其中：λn為第n次倍周期分岔參數(shù)；T0為激勵周期，第n次倍周期分岔后周期Tn=2nT0；x(n)為第n次倍周期分岔后響應(yīng)，n較大時(shí)，周期軌道可近似分成2n個(gè)周期為T0的圈，相鄰兩次倍周期分岔后軌道近似重合。

令：

c(n)(t)=x(n)(t)-x(n)(t+Tn-1)

(2)

式中：上標(biāo)表示第n次倍周期分岔，第n-1次倍周期分岔后周期Tn-1=Tn/2。據(jù)文獻(xiàn)[1-4]對分岔間距收斂常數(shù)α定義[15]：

(3)

據(jù)傅里葉元素定義：

(4)

(5)

聯(lián)合式(3)、(5)推出周期解傅里葉元素奇數(shù)項(xiàng)滿足：

(6)

令2k+1=ξ，2k′+1=ξ′，寫成希爾伯特變換式：

(7)

取P=1，經(jīng)希爾伯特變換后得：

(8)

其中：α=2.502 907 875…為普適分岔間距比，故：

式中：10log10μ=-8.18 dB表示新增線譜平均強(qiáng)度較上次新增線譜平均強(qiáng)度下降8.18 dB。

2 用增量諧波平衡法求Duffing方程周期解

考慮硬彈簧Duffing振子為典型非線性振動模型，其動力學(xué)方程經(jīng)無量綱化后為多項(xiàng)式非線性方程，稱Duffing方程：

(9)

為便于計(jì)算，式(9)改寫為：

(10)

(11)

將式(11)代入式(10)經(jīng)Taylor展開并略去增量高階項(xiàng)得：

(12)

設(shè)周期解初始解x0及增量Δx為：

(13)

(14)

令：

a0=(a0,a1,a2,…,b1,b2,…bNm)

(15)

(16)

Δa0=(Δa0,Δa1,Δa2,…Δb1,Δb2,…ΔbNm)

(17)

則有：

x0=sa0,Δx=sΔa0

(18)

對式(12)應(yīng)用Galerkin過程：

ΔAcos(T+φ)+R]δ(Δx)dT

(19)

其中：

(20)

令：

(21)

(22)

(23)

(24)

(25)

由式(13)～式(25)得：

KΔa0=PΔω+BΔγ+QΔA+R′

(26)

式(26)為 2Nm+1維線性方程組，選Δω，Δγ，ΔA中之一為增量，x0為初始解，除增量外其余各參數(shù)固定；選適當(dāng)增量步長迭代計(jì)算，|R|足夠小時(shí)可求出Δa0，從而求得周期解x(T)。

3 計(jì)算結(jié)果

(27)

第k+1次倍周期分岔,周期2k+1軌道預(yù)表達(dá)式為：

(28)

(29)

(30)

(31)

將式(27)、(28)代替式(13)作為周期解的預(yù)表達(dá)式，經(jīng)上節(jié)計(jì)算過程可求得歷次倍周期分岔后分岔點(diǎn)附近的周期軌道x(k)(T)。周期軌道各譜線幅值可按式(29)計(jì)算，其中Pi(k)為第k次倍周期分岔后第i條譜線幅值強(qiáng)度；由于相鄰兩次倍周期分岔后新增譜線均為奇數(shù)譜線，相鄰兩次倍周期分岔后新增譜線平均強(qiáng)度之比μ按式(30)計(jì)算；分岔參數(shù)收縮率δ按式(31)計(jì)算，其中A(k)為第k次倍周期分岔點(diǎn)激勵幅值。計(jì)算結(jié)果見表1。由表1看出，開始兩次倍周期分岔因不滿足k較大的近似條件，與推導(dǎo)值[1-4]偏差較大，第3、4次倍周期分岔后新增譜線平均強(qiáng)度之比μ4=0.1505已較接近推導(dǎo)值0.1525，且較吻合。

表1 相鄰兩次倍周期分岔后新增譜線平均強(qiáng)度之比μ及分岔參數(shù)收縮率δ

圖1 Duffing振子倍周期分岔的譜特性

將1～6次倍周期分岔后分岔點(diǎn)附近周期解歸一頻率小于1的譜線整合見圖1，圖中1～6標(biāo)號與表1中k對應(yīng)，為1～6次倍周期分岔。第1次倍周期分岔在歸一頻率1/2處增加的譜線為第1條譜線；第2次倍周期分岔在歸一頻率1/4、3/4處增加的兩條譜線為第1、3條譜線，第1次倍周期分岔后兩條譜線仍在原頻率位置；第3次倍周期分岔在1/8、3/8、5/8、7/8處增加四條譜線分別為第1、3、5、7條譜線，上次分岔譜線仍在原頻率位置。每個(gè)倍周期分岔一次，原T倍周期解所有譜線將成為2T倍周期解的偶數(shù)項(xiàng)譜線，新增譜線均為奇數(shù)項(xiàng)譜線。后續(xù)倍周期分岔過程中，偶數(shù)項(xiàng)譜線將留在原頻率位置。式(29)中Pi(k)為第k次倍周期分岔后第i條譜線幅值強(qiáng)度，圖1中橫虛線為第k次倍周期分岔后新增譜線幅值平均強(qiáng)度。第k次倍周期分岔后新增譜線幅值平均強(qiáng)度為：

(32)

(33)

(34)

第k+1次倍周期分岔后新增譜線幅值平均強(qiáng)度較第k次約下降8.2 dB。k較大時(shí)結(jié)果與理論值一致，如第5、6次倍周期分岔后新增譜線幅值平均強(qiáng)度之比μ5取對數(shù)運(yùn)算后結(jié)果為10log10μ6=10log100.152 0=-8.18，亦與實(shí)驗(yàn)觀察結(jié)果吻合。μ6=0.152 0見表1。

4 結(jié) 論

(1)實(shí)驗(yàn)方法只能觀察到前4次分岔，本文用增量諧波平衡法能研究前6次分岔。

(2)與數(shù)值方法相比，增量諧波平衡法所得近似解析解為諧波組合解，無需考慮頻率分辨率，只選傅里葉系數(shù)計(jì)算相鄰兩次倍周期分岔新增譜線平均強(qiáng)度之比μ即可，結(jié)果與理論結(jié)果吻合較好。

(3)本文方法既能驗(yàn)證理論推導(dǎo)的正確性，亦為對數(shù)值、實(shí)驗(yàn)方法驗(yàn)證理論結(jié)論的補(bǔ)充。

參 考 文 獻(xiàn)

[1]Feigenbaum M J. Thetransition to aperiodic behavior in turbulent systems[J].Communications in Mathematical Physics, 1980,77(1): 65-86.

[2]Feigenbaum M J.The onset spectrum of turbulence[J]. Physics Letters A, 1979,74(6):375-378.

[3]Feigenbaum M J. Quantitativeuniversality for a class of nonlinear transformations[J]. Journal of Statistical Physics, 1978,19(1):25-52.

[4]Feigenbaum M J. Universal behavior in nonlinear systems[J]. Los Alamos Science,1980,1:2-27.

[5]Wang L Q,Xu M T.Property of period-doubling bifurcations [J]. Chaos,Solitons and Fractals,2005,24(2):527-532.

[6]Linsay P S. Period doubling and chaotic behavior in a driven anharmonic oscillator[J]. Physical Review Letters,1981, 47(19):1349-1352.

[7]王本仁，繆國慶，魏榮爵.在液氮及水中聲次諧波的實(shí)驗(yàn)觀察[J].物理學(xué)報(bào),1984,33(3):434-436.

WANG Ben-ren,MIAO Guo-qing,WEI Rong-jue.Experimental observation of subharmonic in liquid nitrogen and water[J]. Journal of Acta Physica Sinica, 1984,33(3):434-436.

[8]Tongue B H. Characteristics of numerical simulations of chaotic systems[J]. ASME，Journal of Applied Mechanics，1987,54(3): 695-699.

[9]Skufca J D. Analysis still matters: a surprising instance of failure of Runge-Kutta-Felberg ODEsolvers[J].SIAM Review, 2004,46(4):729-737.

[10]武際可,周 琨.高維極限環(huán)的數(shù)值追蹤[J].中國科學(xué)(A輯)，1994,24(3):269-276.

WU Ji-ke,ZHOU Kun.Numerical tracking of high dimensional limit cycle[J]. Scien in China:Serics A, 1994,24(3):269-276.

[11]Lau S L. The incremental harmonic balance method and its application to nonlinear vibrations[C]. Proceeding of International Conference on Structure Dynamics，Vibration，Noise and Control,Hong Kong,1995:50-57.

[12]Lau S L，Cheung Y K，Wu S Y. Incremental harmonics balance method with multiple time scales for nonlinear aperiodic vibrations[J]. ASME Journal of Applied Mechanics,1983,50: 871-876.

[13]Cheung Y K,Chen S H,Lau S L. Application of the incremental harmonic balance method to cubic nonlinearity systems[J]. Journal of Sound and Vibration,1990,140(2):273-286.

[14]Shen J H,Lin K C,Chen S H, et al. Bifurcation and route-to-chaos analyses for Mathieu-Duffing oscillator by the incremental harmonic balance method[J]. Nonlinear Dynamics,2008,52(4):403-414.

[15]Feigenbaum M J.Presentation functions and scaling function theory for circle maps [J]. Nonlinearity,1988,1(4):577-602.