Duffing振子倍周期分岔譜特性

唐元璋， 樓京俊， 翁雪濤, 朱石堅

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

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

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

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

(1)

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

令：

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

(2)

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

(3)

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

(4)

(5)

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

(6)

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

(7)

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

(8)

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

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

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

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

(9)

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

(10)

(11)

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

(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)增量步長迭代計算，|R|足夠小時可求出Δa0，從而求得周期解x(T)。

3 計算結(jié)果

(27)

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

(28)

(29)

(30)

(31)

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

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

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

將1～6次倍周期分岔后分岔點附近周期解歸一頻率小于1的譜線整合見圖1，圖中1～6標號與表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條譜線，上次分岔譜線仍在原頻率位置。每個倍周期分岔一次，原T倍周期解所有譜線將成為2T倍周期解的偶數(shù)項譜線，新增譜線均為奇數(shù)項譜線。后續(xù)倍周期分岔過程中，偶數(shù)項譜線將留在原頻率位置。式(29)中Pi(k)為第k次倍周期分岔后第i條譜線幅值強度，圖1中橫虛線為第k次倍周期分岔后新增譜線幅值平均強度。第k次倍周期分岔后新增譜線幅值平均強度為：

(32)

(33)

(34)

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

4 結(jié) 論

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

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

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

參 考 文 獻

[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]王本仁，繆國慶，魏榮爵.在液氮及水中聲次諧波的實驗觀察[J].物理學(xué)報,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.