分數階Ver Del Pol-Duffing 系統的非線性動力學行為分析

李媛萍，張 衛(wèi)

(1.重大工程災害與控制教育部重點實驗室(暨南大學)，廣州 510632;2.暨南大學 力學與土木工程系，廣州 510632;3.暨南大學 信息科學技術學院電子工程系，廣州 510632)

經典的Duffing及van der Pol振子雖然在表達形式上很簡單，但是由于具有豐富的動力學特性而極具代表性，它們常常被用來模擬系統的非線性特性。分數導數本構模型由于具有良好的記憶功能，并能在較大頻率范圍內描述材料的力學性能而成為備受重視的本構模型［1，2]，近年來，用分數導數本構模型描述的分數階Duffing振子及分數階van der Pol振子的力學行為已受到科學家們的廣泛關注［3－5]，但是關于分數階 Van der Pol-Duffing耦合振子的動力學行為的研究目前還未見報道?；诖嗽?，本文構造了一類由分數導數本構模型描述的同時含有Van der Pol系統維持自激振動的非線性阻尼項以及Duffing系統的三次非線性恢復力項的分數階Van der Pol-Duffing振子，推導了求解分數導數的數值計算方法，結合Runge-Kutta法對非線性振子方程進行數值求解，探討了該分數階振子在不同類型荷載作用下的動力學特性。

分數階Van der Pol-Duffing振子的狀態(tài)空間描述為:

其中，ε，k＞0為非線性控制參數;ω0為自由系統的固有頻率;F(t)為外部激勵分別為系統任意時刻的位移、速度和加速度;分數導數階值0＜q＜1。當q=1時，系統(1)即為經典的Van der Pol-Duffing系統。

1 數值計算方法

根據Riemann-Liouville分數導數的定義可知，函數f(t)的 q階分數導數為［6]:

設Riemann-Liouville分數導數在tn=n·Δt時刻的值為 Dqf(t)，取等距積分步長 Δt，則 tn=n·Δt，由于被積函數在積分上限τ=tn=n·Δt時有奇性，導致式(2)不能直接求積。當tn足夠大時，如果我們可以把積分區(qū)間分成兩部分:

由于ΔS在積分上限處被積函數奇異，需要進行特殊處理，這里采用等距高斯積分法進行線性化處理，設權重函數為，則:

至此分數導數Dqf(t)即可求出。結合Runge-Kutta法即可對方程(1)進行數值求解。

2 自治系統的動力學特性

首先分析自治系統，即F(t)=0。取初值x1(0)=0.3，x2(0)=0.2，并固定系數 ε =5，ω0=1，k=1，改變分數導數階值q，可以得到不同分數導數階值下系統的時程曲線、功率譜圖及自振周期，分別如圖1和圖2所示。從中可以看出:隨著分數導數階值的變化，系統的振動周期明顯改變，但振幅都趨于極限值2，說明該分數階振子與經典Van der pol振子一樣有極限環(huán)存在也伴有自激瞬態(tài)過程，無論分數導數階值如何改變，當t→∞時，振動將趨于定常振動;從功率譜圖可以看出:該自治系統存在多個振動頻率，假如系統振動的基頻以(0＜q≤1)表示，那么在基頻的奇整數倍kωq1(0＜q≤1，k=3，5，7，…)處將出現多個譜峰，這說明系統的非線性特性非常明顯;同時，隨著分數導數階值的變化，這些譜峰的個數和峰值大小也隨之改變，反映了系統的非線性強弱與分數導數階值密切相關。結合圖2可以看出:系統的振動周期隨參數q的改變而改變，但并非單調遞增或單調遞減的關系，當q=1時系統的自振周期最大，而在分數階區(qū)間，自振周期的變化范圍很大，其中q=0.4時系統的自振周期最小，僅為q=1時的45%，在q∈(0.3～0.5)區(qū)間，系統的自振周期比較接近，都比q=1時的最大值小了50%以上，由此可見，通過改變分數導數階值的大小可以調節(jié)系統的自振特性，而且調節(jié)的范圍還很大，這一特點正體現了分數導數本構模型的優(yōu)越性。由于分數導數階值的變化實質上反映了系統阻尼特性的改變，從上面的分析可知，在參數q∈(0.3～0.5)范圍內系統具有較好的阻尼性能。

圖1 q=1，0.5，0.1 時自治系統的功率譜Fig.1 Power spectra for autonomous system with fractional order q={1，0.5，0.1}

圖 2 q=0.1∶0.1∶1 時自治系統的振動周期Fig.2 period for autonomous system with fractional order q=0.1∶0.1∶1

維持初值 x1(0)=0.3，x2(0)=0.2 及系數 ε =1，ω0=1，q=0.5不變，改變非線性剛度系數 k，經數值計算可得到不同非線性剛度系數k下自治系統的相圖、時域曲線和功率譜，見圖3。從圖中可以看出，隨著非線性剛度系數的增大，極限環(huán)的形狀和大小均發(fā)生了明顯改變，同時系統的振動周期逐漸減小，但振幅仍限定在極限值2以內，說明無論非線性剛度系數k如何增大，系統始終能趨于定常振動且該定常振動是非常穩(wěn)定的;同時，從功率譜圖可以明顯看到，隨著非線性剛度系數k的增大，系統除在基頻和少數奇數倍頻處有尖峰外，還出現了分頻成分，而且k值越大分頻成分越多，反映了系統的非線性特性隨k值的增大而顯著增強。

圖3 k=0，5，10時自治系統的相圖、時域曲線和功率譜Fig.3 Phase diagram(left)and time domain(middle)and power spectra(right)with cubic nonlinear resilience coefficient k={0，5，10}

圖4 ε=1.5，10時自治系統的相圖、時間位移曲線和功率譜Fig.4 Phase diagram(left)and time domain(middle)and power spectra(right)with damping coefficientε =1.5，10

維持初值 x1(0)=0.3，x2(0)=0.2 及系數 ω0=1，q=0.5，k=1不變，改變阻尼系數ε，可得到系統自振特性隨阻尼系數的變化規(guī)律，如圖4所示。從圖中可以看出:隨著阻尼系數的增大，極限環(huán)圍繞的面積增大且形狀逐漸變得歪扭，功率譜曲線上顯示系統振動的倍頻分量逐漸增多，同時時間位移曲線的波形逐漸由簡諧振動變?yōu)閺埑谡駝?，但振幅始終限定在極限值2以內，反映系統的非線性特性隨著阻尼系數的增大而增強，但無論阻尼如何變化，當t→∞時，系統都將趨于定常振動，這些正是該分數階振子與經典Van der Pol-Duffing振子的相似之處。

3 周期載荷下的動力學特性

考慮外部激勵為簡諧荷載F(t)=p cos(ωt)的情形。取參數 ω =1.2，ε =0.1，q=0.5，k1=1=1，同時取系統的初始條件為x1(0)=x2(0)=0，經數值計算可得到系統的位移隨外激勵幅值變化的分岔圖，如圖5所示。從圖中可以看出:隨著外激勵幅值的增大，系統將由擬周期振動(見圖6)變?yōu)橹芷谌駝?見 圖7)，最后發(fā)展為單周期振動(見圖8)。在該范圍內未發(fā)現混沌現象。

圖5 系統的位移隨參數p變化的分岔圖Fig.5 bifurcation diagram of displacement with varying of parameter p at scope of(1.5，2.5)

圖6 外激勵幅值p=1.6時系統的相軌跡和Poincare截面Fig.6 Phase portrait and Poincare map at p=1.6

圖7 外激勵幅值p=2.3時系統的相軌跡和Poincare截面Fig.7 Phase portrait and Poincare map at p=2.3

圖8 外激勵幅值p=2.5系統的相軌跡和Poincare截面Fig.8 Phase portrait and Poincare map at p=2.5

當外部荷載為簡諧荷載時，在系統中盡管存在著非線性的影響，但在一定條件下其定常解仍近似于簡諧的，基于此，我們可以利用傅里葉變換對其諧波成分進行定量分析。圖9是在外激勵幅值p∈(1.5，2.5)區(qū)間內各次諧波分量的幅值隨外激勵幅值的變化曲線，計算中共取了五次諧波加以分析，其中0次諧波表示a0的值，從圖中我們可以清楚地看出:在擬周期振動區(qū)間，前三次諧波為主要成分且幅值波動很大;在單周期振動區(qū)間，主要成分為一次諧波;在周期三振動區(qū)間，前三次諧波分量均占有較大的比例且波動較小。從圖中還可以看出，在整個變化區(qū)間，第四、五次諧波分量的幅值均比較小，說明其所占的比例很小，為了簡化計算，我們完全可以只取前三次諧波加以分析。

圖9 各諧波分量的幅值隨外激勵幅值的變化曲線Fig.9 Change of harmonic component with excitation amplitude

圖10 系統的位移隨阻尼系數ε變化的分岔圖Fig.10 Bifurcation diagram of displacement with varying of damping coefficient at scope ofε

下面探討在簡諧激勵下阻尼系數的變化對系統振動特性的影響。假設外激勵的幅值和頻率分別為p=2.5，ω =1.2，同時假設系統的初始狀態(tài)為 x1(0)=x2(0)=0，固定參數 q=0.5，k1=1=1，可得到系統位移隨阻尼系數ε變化的分岔圖，如圖10所示。從圖中可見，隨著阻尼系數的增大，系統將由單周期振動(圖11)變?yōu)橹芷谌駝?圖12)最后發(fā)展擬周期振動(圖13)。在該范圍內未發(fā)現混沌現象。

圖11 阻尼系數ε=0.15時系統的相軌跡和Poincare截面Fig.11 Phase portrait and Poincare map atε =0.15

圖12 阻尼系數ε=0.95時系統的相軌跡和Poincare截面Fig.12 Phase portrait and Poincare map atε =0.95

圖13 阻尼系數ε=1時系統的相軌跡和Poincare截面Fig.13 Phase portrait and Poincare map atε =1

4 地震荷載下的運動特性

最后分析系統在地震荷載作用下的振動特性。地震波的實質是一種非平穩(wěn)信號，因此在地震作用下結構的反應無疑也是非平穩(wěn)的，從而使得在非線性分析中普遍采用的相圖、Poincare截面以及分岔圖等分析方法已不再適用，本文采用信號處理中常用的快速傅里葉變換法進行處理。選取EL－centro地震波的南北分量作為輸入(加速度峰值為341.7 m/s2，取樣時間為前30 s，取樣頻率為50 Hz)，其加速度時程曲線和相應的功率譜見圖14，從其功率譜圖可看出，其能量主要分布在0 Hz～5 Hz的范圍內。

固定參數 ε =0.1，k1=1=1，并設系統的初始狀態(tài)為x1(0)=x2(0)=0，改變q，可得到系統在不同分數導數階值下的動力反應特征。為節(jié)省篇幅，本文僅示意出分數導數階值分別為 q=0.1，0.5，0.9 時系統動力反應的時域曲線和功率譜，見圖15。從圖中可以看出:EL－centro南北波經過該系統之后，輸出信號的頻率范圍和幅值均大大減小且主要集中于0.5 Hz～1 Hz范圍內，說明大部分的輸入能量被消耗;其中，當q=0.5時系統輸出信號的功率譜的幅值最小，同時其位移幅值隨著時間的增加逐漸衰減，說明q=0.5時系統輸出的動力反應最小，究其原因正是因為當q=0.5時系統的阻尼較大，從而消耗掉更多的能量。

圖14 EL-centro地震波時程曲線及功率譜Fig.14 Time domain and power spectra of EL-centro wave

圖15 q=0.1，0.3，0.5，0.7，0.9 時系統動力反應的時域曲線和功率譜Fig.15 Time domain and power spectra with fractional order q=0.1，0.3，0.5，0.7，0.9 under seismic load

5 結論

本文引進分數導數本構模型模擬系統的阻尼特性，構造了一類同時含有Van der Pol系統維持自激振動的非線性阻尼項以及Duffing系統的三次非線性恢復力項的分數階Van der Pol-Duffing振子，通過數值計算探討了該系統在不同荷載作用下的振動特性，可以得到如下結論:

(1)該非線性振子具有與經典Ver Del Pol系統相似的自激振動特性，同時分數導數階值的變化能改變系統的自振周期;在q∈(0.3～0.5)區(qū)間，系統的自振周期比整數階系統減小一半以上。

(2)隨著阻尼系數和非線性位移系數的增大，系統的非線性增強。

(3)在簡諧荷載作用下，隨著外荷載幅值的增大，系統由擬周期振動變?yōu)橹芷谌駝幼詈蟀l(fā)展為單周期振動。

(4)在簡諧荷載作用下，隨著阻尼系數的增大，系統由單周期振動變?yōu)橹芷谌駝幼詈蟀l(fā)展為擬周期振動。

(5)在地震荷載作用下，通過調節(jié)分數導數階值可以改變系統的阻尼性能從而改變輸出的能量。

［1] Ge Z M，Ou C Y，Chaos in a fractional order modified duffing system［J] .Chaos Solitons ＆ Fractals，2007，34(2):262－291.

［2] 李根國，朱正佑，程昌鈞.非線性粘彈性Timoshenko梁動力學行為的分析［J] .力學季刊，2001.22(3):346－352.

［3] Gao X，Yu J B.Chaos in the fractional order periodically forced complex Duffing's oscillators［J] .Chaos Solitons ＆Fractals，2005.24(4):1097 －1104.

［4] Cao J，Xie H，Jiang Z.Nonlinear dynamics of duffing system with fractional order damping［J] .Hsi-An Chiao Tung Ta Hsueh/Journal of Xi'an Jiaotong University，2009，43(3):50－54.

［5] Barbosa R S，Machado J A T，Ferreira I M.Dynamics of the fractional-order van der pol Oscillator［C] .IEEE Transactions on Industrial Electronics:373－378.

［6] 張 衛(wèi)，徐 華，清水信行.分數算子描述的粘彈性體力學問題數值方法［J] .力學學報，2004，36(5):617－622.