李媛萍，張 衛(wèi)

(1.重大工程災(zāi)害與控制教育部重點(diǎn)實(shí)驗(yàn)室(暨南大學(xué))，廣州 510632;2.暨南大學(xué) 力學(xué)與土木工程系，廣州 510632;3.暨南大學(xué) 信息科學(xué)技術(shù)學(xué)院電子工程系，廣州 510632)

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

分?jǐn)?shù)階Van der Pol-Duffing振子的狀態(tài)空間描述為:

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

1 數(shù)值計算方法

根據(jù)Riemann-Liouville分?jǐn)?shù)導(dǎo)數(shù)的定義可知，函數(shù)f(t)的 q階分?jǐn)?shù)導(dǎo)數(shù)為［6]:

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

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

至此分?jǐn)?shù)導(dǎo)數(shù)Dqf(t)即可求出。結(jié)合Runge-Kutta法即可對方程(1)進(jìn)行數(shù)值求解。

2 自治系統(tǒng)的動力學(xué)特性

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

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

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

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

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

3 周期載荷下的動力學(xué)特性

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

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

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

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

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

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

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

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

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

圖11 阻尼系數(shù)ε=0.15時系統(tǒng)的相軌跡和Poincare截面Fig.11 Phase portrait and Poincare map atε =0.15

圖12 阻尼系數(shù)ε=0.95時系統(tǒng)的相軌跡和Poincare截面Fig.12 Phase portrait and Poincare map atε =0.95

圖13 阻尼系數(shù)ε=1時系統(tǒng)的相軌跡和Poincare截面Fig.13 Phase portrait and Poincare map atε =1

4 地震荷載下的運(yùn)動特性

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

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

圖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 時系統(tǒng)動力反應(yīng)的時域曲線和功率譜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 結(jié)論

本文引進(jìn)分?jǐn)?shù)導(dǎo)數(shù)本構(gòu)模型模擬系統(tǒng)的阻尼特性，構(gòu)造了一類同時含有Van der Pol系統(tǒng)維持自激振動的非線性阻尼項(xiàng)以及Duffing系統(tǒng)的三次非線性恢復(fù)力項(xiàng)的分?jǐn)?shù)階Van der Pol-Duffing振子，通過數(shù)值計算探討了該系統(tǒng)在不同荷載作用下的振動特性，可以得到如下結(jié)論:

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

(2)隨著阻尼系數(shù)和非線性位移系數(shù)的增大，系統(tǒng)的非線性增強(qiáng)。

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

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

(5)在地震荷載作用下，通過調(diào)節(jié)分?jǐn)?shù)導(dǎo)數(shù)階值可以改變系統(tǒng)的阻尼性能從而改變輸出的能量。

［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梁動力學(xué)行為的分析［J] .力學(xué)季刊，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，F(xiàn)erreira I M.Dynamics of the fractional-order van der pol Oscillator［C] .IEEE Transactions on Industrial Electronics:373－378.

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