分數階擬周期Mathieu 方程的動力學分析1)

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

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

? (石家莊鐵道大學交通工程結構力學行為與系統安全國家重點實驗室,石家莊 050043)

引言

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

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

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

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

1 分數階擬周期Mathieu 方程

1.1 攝動分析

本文研究的分數階擬周期Mathieu 方程如下所示

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

關于分數階微積分的定義有多種,這里采用Caputo 型分數階微積分定義

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

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

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

過渡曲線在 δ-ω 平面內具有如下形式

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

由式(6a)解得

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

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

根據式(6a)解得

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

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

則式(9)的特解為

這里利用公式[30]對分數階微分項kDpu1進行處理

再結合歐拉公式得到

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

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

此時式(14)的特解為

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

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

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

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

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

對分數階項kDpu0進行處理得到

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

消除永年項條件為

該方程有非零解的條件是

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

此時式(6b)的特解為

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

為了消除永年項,要求

從而得到

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

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

1.2 方程過渡曲線分析

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

圖1 分數階擬周期Mathieu 方程的過渡曲線Fig.1 Transition curves of QP Mathieu equation with fractional-order derivative

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

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

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

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

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

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

2 數值仿真

2.1 數值解和解析解的比較

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

其中tl=lh為時間采樣點,h為時間步長,為分數階二項式系數,具有以下迭代關系

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

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

2.2 分數階微分項對過渡曲線的影響

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

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

圖3 分數階微分項階次 p 對過渡曲線的影響Fig.3 Effects of the fractional order p on transition curves

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

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

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

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

圖5 p=0.5 時分數階微分項系數對過渡曲線的影響Fig.5 Effects of the fractional coefficient on transition curves when p=0.5

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

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

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

3 結論

應用攝動法研究了分數階擬周期Mathieu 方程,得到了方程在 δ-ω 平面內過渡曲線的近似表達式.借助等效線性剛度和等效線性阻尼概念,分別分析了不同分數階微分項系數和階次對方程過渡曲線的影響.結果發現,分數階微分項同時具有剛度特性和阻尼特性,選取不同的分數階微分項階次和系數可以使其呈現不同程度的剛度特性或阻尼特性,方程穩定區域的大小和過渡曲線的位置也因此產生了不同程度的變化.以上結果說明分數階微分項對擬周期Mathieu 方程的穩定特性有著重要的影響,對此類系統的分析和穩定狀態參數的選取有著重要的意義.