郭建斌 申永軍 ,?,2) 李 航

* (石家莊鐵道大學(xué)機械工程學(xué)院,石家莊 050043)

? (石家莊鐵道大學(xué)交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國家重點實驗室,石家莊 050043)

引言

分?jǐn)?shù)階微積分第一次被數(shù)學(xué)家Hospital 與Leibnitz 提出,至今已有300 多年歷史.起初,分?jǐn)?shù)階微積分因其缺乏應(yīng)用背景在早期工程研究中并未得到廣泛關(guān)注.1832 年,Liouville 首次將分?jǐn)?shù)階微積分用于解決一些實際問題.此后,眾多學(xué)者在分?jǐn)?shù)階微積分的定義、性質(zhì)和計算方法等方面展開了詳細(xì)的探討[1-6].在這個過程中,分?jǐn)?shù)階微積分也由數(shù)學(xué)理論研究逐步走向工程應(yīng)用[7-11].

在控制工程領(lǐng)域,分?jǐn)?shù)階微積分主要用來影響系統(tǒng)的閉環(huán)控制特性,從而提高系統(tǒng)控制效果和魯棒性.如薛定宇和趙春娜[12]提出了一種分?jǐn)?shù)階PID控制器的設(shè)計方法,具體演示了分?jǐn)?shù)階控制器對系統(tǒng)優(yōu)良的控制效果.常宇健等[13]則針對含分?jǐn)?shù)階阻尼的懸架模型進(jìn)行了主動控制研究.在動力學(xué)方面,分?jǐn)?shù)階微積分主要集中用于描述黏彈特性材料的本構(gòu)關(guān)系,以此來提高此類非線性系統(tǒng)振動特性研究的準(zhǔn)確性.如Cao 等[14]采用數(shù)值積分法,結(jié)合相圖、龐加萊截面圖、分岔圖等分析了分?jǐn)?shù)階阻尼對系統(tǒng)動力學(xué)性能的影響.申永軍等[15]通過對含分?jǐn)?shù)階線性單自由度振子的動力學(xué)分析,首次提出等效線性阻尼和等效線性剛度概念.申永軍等[16-17]利用平均法得到分?jǐn)?shù)階van der Pol 振子的一次近似解,比較了分?jǐn)?shù)階與整數(shù)階系統(tǒng)并分析了分?jǐn)?shù)階微分項對系統(tǒng)動力特性的影響.由上述研究可知,相比整數(shù)階模型,分?jǐn)?shù)階模型的物理意義更清晰,更能準(zhǔn)確描述實際系統(tǒng).

Mathieu 方程作為一種周期系數(shù)線性微分方程,因其本身復(fù)雜的動力學(xué)特性,在動力學(xué)領(lǐng)域得到了廣泛的應(yīng)用[18-20],經(jīng)常被用來處理一些參數(shù)共振現(xiàn)象.如Qian 等[21]研究了在參數(shù)激勵和強迫激勵作用下的斜拉橋拉索的非線性動力學(xué)問題,利用多尺度方法分析了1/2 參激共振下的動力學(xué)響應(yīng).黃建亮等[22]對van der Pol-Mathieu 方程的動力學(xué)特性進(jìn)行了研究,應(yīng)用改進(jìn)的諧波平衡法精確計算出了方程的準(zhǔn)周期響應(yīng).溫少芳[23]研究了分?jǐn)?shù)階參激系統(tǒng)的動力學(xué)特性,分析了分?jǐn)?shù)階微分項對系統(tǒng)穩(wěn)定性邊界以及幅頻響應(yīng)的影響.同時,隨著眾多學(xué)者對Mathieu 方程的深入研究,促使其形式也得到了豐富拓展.Kovacic 等[24]對Mathieu 方程的推廣型作了全面概述,包含擬周期系數(shù)和橢圓系數(shù)等多種形式的Mathieu 方程.其中,擬周期系數(shù)Mathieu 方程通常被應(yīng)用于一些特殊的動力系統(tǒng),如Galeotti 和Toni[25]提出的含雙頻激勵時變剛度的高速列車受電弓非線性模型.Huan 等[26]研究了高速列車受電弓-懸鏈線組合系統(tǒng)中低頻和高頻參數(shù)激勵對系統(tǒng)響應(yīng)的影響等.Rand 等[27-28]則針對擬周期Mathieu 方程不同共振狀態(tài)的穩(wěn)定區(qū)域進(jìn)行了深入研究.

綜上所述,在含有黏彈性器件的參激系統(tǒng)(如弓網(wǎng)系統(tǒng))中,分?jǐn)?shù)階模型較整數(shù)階模型對系統(tǒng)描述更加準(zhǔn)確.因此,本文在擬周期Mathieu 方程中引入分?jǐn)?shù)階微積分理論,建立了分?jǐn)?shù)階擬周期Mathieu 方程,利用攝動法求得了方程過渡曲線的近似解析解,通過數(shù)值方法在系統(tǒng)穩(wěn)定圖中分析了分?jǐn)?shù)階微分項參數(shù)對方程穩(wěn)定區(qū)域大小和過渡曲線位置的影響,并驗證了分?jǐn)?shù)階微分項同時含有阻尼和剛度特性的特點.

1 分?jǐn)?shù)階擬周期Mathieu 方程

1.1 攝動分析

本文研究的分?jǐn)?shù)階擬周期Mathieu 方程如下所示

其中,ω 是無理數(shù);ε 為小參數(shù),滿足| ε|?1 ;2 ζ 和δ+ε[cost+cos(ωt)]分別為線性阻尼和時變剛度;Dp[u(t)]為u(t) 關(guān)于t的p階導(dǎo)數(shù) (0 ≤p≤1);K1為分?jǐn)?shù)階微分項系數(shù) (K1＞0) .

關(guān)于分?jǐn)?shù)階微積分的定義有多種,這里采用Caputo 型分?jǐn)?shù)階微積分定義

式中 Γ 為Gamma 函數(shù),具有 Γ(z+1)=zΓ(z) 的特性.

為了確定方程穩(wěn)態(tài)周期解的過渡曲線,引入ζ=εμ,μ=O(1),K1=εk,k=O(1),式(1)變換為

利用攝動法,假設(shè)式(3)的解滿足

過渡曲線在 δ-ω 平面內(nèi)具有如下形式

將式(4)和式(5)代入式(3),比較 ε 的同次冪得到

由式(6a)解得

其中cc代表前面各項的共軛,后不贅述.由于 ω 是無理的,擬周期Mathieu 方程穩(wěn)態(tài)周期解的過渡曲線在δ-ω平面上是從處產(chǎn)生[29],下面討論方程在 α=0,1,β=-1,0,1 時的過渡曲線.

(1) δ0=0 (即 α=0,β=0)

根據(jù)式(6a)解得

式中c是由初始條件確定的常量,將其代入式(6b)得到

為了消除永年項需 -δ1c=0,所以有

則式(9)的特解為

這里利用公式[30]對分?jǐn)?shù)階微分項kDpu1進(jìn)行處理

再結(jié)合歐拉公式得到

將式(8)、式(10)、式(11)和式(13)代入式(6c)得到

從上式得出消除永年項條件

此時式(14)的特解為

將式(10)和式(15)代入式(5)得到

分析式(17)可以看出,在 δ0=0 附近過渡曲線的二階近似解與方程的分?jǐn)?shù)階微分項和阻尼都沒有關(guān)系,所以進(jìn)行三階近似計算

類似地,將式(8)、式(10)、式(11)、式(15)和式(16)代入式(18)可得到消除永年項條件

代入方程原參數(shù)得到 δ0=0 時過渡曲線的三階近似解

由式(7)得到式(6a)的特解為

對分?jǐn)?shù)階項kDpu0進(jìn)行處理得到

將式(21)和式(22)代入式(6b)中得到

消除永年項條件為

該方程有非零解的條件是

其中det 為求解矩陣行列式.計算得出

此時式(6b)的特解為

將式(21)、式(27)和式(28)代入式(6c)得

為了消除永年項,要求

從而得到

將式(27)和式(30)代入式(5),代入原參數(shù)整理得到此時的兩條過渡曲線

根據(jù)上述類似的計算方法相應(yīng)求出β=0)時方程的過渡曲線為

1.2 方程過渡曲線分析

選取一組參數(shù) ζ=0.005,K1=0.005,p=0.5,ε=0.1,利用式(20)和式(31)～ 式(34)并略去式中高階部分繪制解析結(jié)果的過渡曲線如圖1 所示.

圖1 分?jǐn)?shù)階擬周期Mathieu 方程的過渡曲線Fig.1 Transition curves of QP Mathieu equation with fractional-order derivative

對比不含分?jǐn)?shù)階微分項的Mathieu 系統(tǒng)[31-32],通過定義等效線性阻尼C(p) 和等效線性剛度K(p) 的方法來分析分?jǐn)?shù)階微分項對擬周期Mathieu 方程過渡曲線的影響. δ0在各情況下的等效線性阻尼和等效線性剛度如表1 所示.

通過表1,得到方程等效阻尼和等效剛度的一般形式

表1 δ0 在不同情況下的等效線性阻尼和等效線性剛度Table 1 Equivalent linear damping and equivalent linear stiffness of different δ0

分析上述5 種情況下的結(jié)果可知,在 δ0=0 時,方程過渡曲線的二階近似解與分?jǐn)?shù)階微分項無關(guān),分?jǐn)?shù)階微分項以等效線性阻尼和等效線性剛度的形式影響著過渡曲線的三階近似解.而在其他4 種情況中,分?jǐn)?shù)階微分項對方程二階近似解均有影響,并且它們的等效線性阻尼和等效線性剛度均可整理為一般形式(35).通過分析式(35)發(fā)現(xiàn),分?jǐn)?shù)階微分項的系數(shù)K1和階次p對方程過渡曲線有著重要影響:當(dāng)分?jǐn)?shù)階微分項系數(shù)K1逐漸增大時,等效線性阻尼和等效線性剛度都會逐漸增大;當(dāng)分?jǐn)?shù)階階次p趨近于0 時,分?jǐn)?shù)階微分項幾乎等于線性剛度;而當(dāng)p趨近于1 時,分?jǐn)?shù)階微分項幾乎等于線性阻尼.

另外,通過對比式(31)～ 式(34)發(fā)現(xiàn)方程過渡曲線表達(dá)式具有一定的相似性,都是由 δ0開始隨著ω的增大逐漸分裂成兩條過渡曲線組成(δ=δ0+E±F),如方程在時過渡曲線式(31).根據(jù)這一特點,可以定義此時兩條過渡曲線之間的厚度(即非穩(wěn)定區(qū)域厚度)為

類似地,在 δ0為其他幾種情況時也存在類似的式子來表示過渡曲線之間的厚度.利用式(36)可以更加直觀地顯示分?jǐn)?shù)階微分項系數(shù)和階次的變化對δ-ω平面上非穩(wěn)定性區(qū)域大小的影響.

2 數(shù)值仿真

2.1 數(shù)值解和解析解的比較

為了驗證本文結(jié)果的正確性,下面將上述過渡曲線的解析結(jié)果和數(shù)值結(jié)果進(jìn)行對比.利用文獻(xiàn)[5]中介紹的數(shù)值方法研究方程(1),該方法的近似公式為

其中tl=lh為時間采樣點,h為時間步長,為分?jǐn)?shù)階二項式系數(shù),具有以下迭代關(guān)系

將 ω=0～1.5 和δ=-0.2～0.5 組成的δ -ω 平面內(nèi)所有的點離散處理,分別代入式(1)中對方程進(jìn)行數(shù)值積分,計算一段時間后根據(jù)方程響應(yīng)的振幅變化情況判斷各個參數(shù)點對應(yīng)的穩(wěn)定性,以此來確定式(1)的穩(wěn)定區(qū)和非穩(wěn)定區(qū)分界線.其中 δ 選擇間隔為0.002,ω 間隔為0.01,計算時間為700 s,計算步長選擇0.001,所得結(jié)果如圖2 所示.仿真過程中參數(shù)取值為: ζ=0.005,K1=0.005,p=0.5,ε=0.1 .圖中黑色的點代表數(shù)值解穩(wěn)定點,白色是非穩(wěn)定區(qū)域,白色和黑色的分界線是數(shù)值解的穩(wěn)定性邊界,紅圈代表方程過渡曲線的解析結(jié)果.從圖2 中可以看到方程在附近形成了指狀的非穩(wěn)定區(qū)域,且解析結(jié)果和數(shù)值結(jié)果的穩(wěn)定性邊界在主要區(qū)域吻合度較好,證明了文中所述方法和結(jié)果具有較好的準(zhǔn)確性.

圖2 數(shù)值解和解析解的方程過渡曲線Fig.2 Transition curves of numerical and analytical solutions

2.2 分?jǐn)?shù)階微分項對過渡曲線的影響

下面研究分?jǐn)?shù)階微分項參數(shù)對方程過渡曲線的影響.首先選定參數(shù) ε,K1和 ζ,當(dāng)分?jǐn)?shù)階微分項階次p分別選取0.1,0.5 和0.9 時,方程在時的過渡曲線如圖3 所示.

從圖3 中可以看出,當(dāng)分?jǐn)?shù)階微分項階次p逐漸增大時,由于等效線性剛度K(p) 在逐漸減小,因此方程的過渡曲線在向右移動;同時,等效線性阻尼C(p)在逐漸增大,方程的非穩(wěn)定區(qū)域在逐漸縮小.說明分?jǐn)?shù)階階次p不僅影響方程過渡曲線的位置,而且還影響方程穩(wěn)定區(qū)域的大小.

圖3 分?jǐn)?shù)階微分項階次 p 對過渡曲線的影響Fig.3 Effects of the fractional order p on transition curves

為了更好地區(qū)別分?jǐn)?shù)階微分項的阻尼和剛度特性,分別取分?jǐn)?shù)階微分項階次p=0.1,p=0.5和p=0.9,觀察不同階次對方程過渡曲線的影響(ε=0.1,ζ=0.005),所得結(jié)果如圖4～ 圖6 所示.

圖4 中給出了p=0.1 時方程過渡曲線隨系數(shù)K1的變化情況.可以發(fā)現(xiàn),當(dāng)K1逐漸增大時,由于等效線性剛度也在逐漸增大,因此方程過渡曲線的位置發(fā)生了明顯的左偏移;但此時等效線性阻尼變化較小,所以非穩(wěn)定區(qū)域的面積未出現(xiàn)明顯收縮.說明p=0.1時分?jǐn)?shù)階微分項呈現(xiàn)出較強的剛度特性,而阻尼特性相對較弱.

圖4 p=0.1時分?jǐn)?shù)階微分項系數(shù) K1 對過渡曲線的影響Fig.4 Effects of the fractional coefficient K1 on transition curves when p=0.1

當(dāng)分?jǐn)?shù)階微分項階次p=0.5 時,如圖5(a) 和圖5(b)顯示,隨著分?jǐn)?shù)階微分項系數(shù)K1的逐漸增大,此時等效線性剛度和等效線性阻尼都在增大,方程過渡曲線不僅逐漸向左偏移,同時曲線之間的厚度也在逐漸縮小.利用式(36) 觀察時非穩(wěn)定區(qū)厚度隨系數(shù)K1的變化情況,從圖5(c)看出,隨著系數(shù)K1的逐漸增大,方程非穩(wěn)定區(qū)面積在逐漸縮小,并且當(dāng)系數(shù)K1增大到一定程度,發(fā)生了局部非穩(wěn)定區(qū)域消失的現(xiàn)象.說明當(dāng)p=0.5 時,分?jǐn)?shù)階微分項不僅具有明顯的剛度特性還具有明顯的阻尼特性.

圖5 p=0.5 時分?jǐn)?shù)階微分項系數(shù)對過渡曲線的影響Fig.5 Effects of the fractional coefficient on transition curves when p=0.5

圖6 p=0.9 時分?jǐn)?shù)階微分項系數(shù) K1 對過渡曲線的影響Fig.6 Effects of the fractional coef0ficient K1 on transition curves when p=0.9

在圖6 中分析p=0.9 時不同分?jǐn)?shù)階微分項系數(shù)K1對過渡曲線的影響.此外,為了便于比較,給出了p=0.5,K1=0.001時線性阻尼系數(shù) ζ 對方程過渡曲線的影響情況,如圖7 所示.由圖6 和圖7 可知,隨著系數(shù)K1和 ζ 逐漸增大,方程的非穩(wěn)定區(qū)域都發(fā)生了明顯的收縮,說明此時分?jǐn)?shù)階微分項呈現(xiàn)出較強的阻尼特性,系數(shù)K1對系統(tǒng)的作用與線性阻尼系數(shù) ζ 幾乎相同.

圖7 線性阻尼系數(shù) ζ 對過渡曲線的影響Fig.7 The evolutions of the transition curves due to the change ofζ

3 結(jié)論

應(yīng)用攝動法研究了分?jǐn)?shù)階擬周期Mathieu 方程,得到了方程在 δ-ω 平面內(nèi)過渡曲線的近似表達(dá)式.借助等效線性剛度和等效線性阻尼概念,分別分析了不同分?jǐn)?shù)階微分項系數(shù)和階次對方程過渡曲線的影響.結(jié)果發(fā)現(xiàn),分?jǐn)?shù)階微分項同時具有剛度特性和阻尼特性,選取不同的分?jǐn)?shù)階微分項階次和系數(shù)可以使其呈現(xiàn)不同程度的剛度特性或阻尼特性,方程穩(wěn)定區(qū)域的大小和過渡曲線的位置也因此產(chǎn)生了不同程度的變化.以上結(jié)果說明分?jǐn)?shù)階微分項對擬周期Mathieu 方程的穩(wěn)定特性有著重要的影響,對此類系統(tǒng)的分析和穩(wěn)定狀態(tài)參數(shù)的選取有著重要的意義.