Duffing振子倍周期分岔譜特性

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

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

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

1 倍周期分岔線譜變化規律

設N維非線性連續時間動力系統為：

(1)

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

令：

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

(2)

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

(3)

據傅里葉元素定義：

(4)

(5)

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

(6)

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

(7)

取P=1，經希爾伯特變換后得：

(8)

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

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

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

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

(9)

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

(10)

(11)

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

(12)

設周期解初始解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)應用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為初始解，除增量外其余各參數固定；選適當增量步長迭代計算，|R|足夠小時可求出Δa0，從而求得周期解x(T)。

3 計算結果

(27)

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

(28)

(29)

(30)

(31)

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

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

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

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

(32)

(33)

(34)

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

4 結 論

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

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

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

參 考 文 獻

[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].物理學報,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]武際可,周 琨.高維極限環的數值追蹤[J].中國科學(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.