分數(shù)階非線性隔振系統(tǒng)的超諧波共振與周期運動轉(zhuǎn)遷規(guī)律分析

屈鳴鶴， 吳少培， 俞力洋， 丁旺才， 李國芳， 黃 然

(蘭州交通大學 機電工程學院，蘭州 730070)

黏彈性隔振元件被廣泛應用于機械設備中，如航空器的輔助動力裝置隔振器[1]、軌道車輛二系懸掛的橡膠堆、汽車底盤中的懸架緩沖塊[2]、機床的減震墊等。其中橡膠只有在變形較小時才可近似看作線性彈性材料，超出這個范圍表現(xiàn)為非線性彈性。由于其特殊的結(jié)構(gòu)可以同時表現(xiàn)出黏彈性，能集緩沖、隔振、降噪等功能于一身，并且可根據(jù)實際需要來設計外形、剛度及阻尼[3]。為了更好地描述黏彈性隔振系統(tǒng)的力學特性，有學者提出三元件固體整數(shù)階模型，可以更好地反映黏彈性材料的松弛和蠕變特性[4]。學者們將可以等效橡膠材料的三參數(shù)力學模型稱為Zener模型和Ruzicka模型[5-7]。

分數(shù)階微積分由于缺乏實際應用背景使其長期以來沒有得到研究和發(fā)展，近幾年分數(shù)階微積分的定義、特性和計算才得以在工程領(lǐng)域應用[8]。雖然傳統(tǒng)的整數(shù)階模型可以描述橡膠材料的力學特性，但不足以描述黏彈材料的頻率相關(guān)性。分數(shù)階力學模型可以對黏彈性加以修正，進而將整數(shù)階模型優(yōu)化為分數(shù)階模型描述材料的本構(gòu)關(guān)系[9]。羅文波等[10]為準確描述瀝青混合料的動態(tài)黏彈性力學行為，在分數(shù)階Zener模型的基礎上提出了改進的分數(shù)階Zener模型。Martin等[11]對黏彈性納米梁的分數(shù)階動力學行為進行研究，借助拉普拉斯變換等方法對黏彈性特性進行了研究分析。Bratu等[12]通過分數(shù)階Zener模型描述了道路的復合土結(jié)構(gòu)，并分析了系統(tǒng)的動態(tài)響應。Lewandowski等[13]用分數(shù)階Zener模型來反應材料的流變特性。Ciniello等[14]分析了溫度對黏彈性材料的影響。

針對隔振系統(tǒng)的非線性特性，常宇健等[15]提出一種含有分數(shù)階微分的金屬橡膠黏彈性本構(gòu)模型，在此基礎上建立了非線性動力學模型，并結(jié)合實驗驗證了模型的準確性。秦浩等對Caputo定義下的分數(shù)階Duffing振子解析解和數(shù)值解的比較，驗證了分數(shù)階項化簡為一階三角函數(shù)形式系統(tǒng)產(chǎn)生分岔和混沌的必要條件。孔凡等[16]通過諧波平衡法研究了簡諧激勵下同時具有滯回特性和分數(shù)階阻尼單元系統(tǒng)的穩(wěn)態(tài)響應，采用不同方法求解系統(tǒng)的遲滯回線，發(fā)現(xiàn)諧波平衡法與逐步積分法得到的相關(guān)結(jié)果吻合較好。Zhen等[17]通過高次諧波平衡法求解了含負剛度幾何非線性系統(tǒng)的動力學響應，并對其隔振性能進行了分析。余慧杰等[18]通過三次非線性函數(shù)描述金屬橡膠的非線性特性，用分數(shù)階模型描述其黏彈性，所建立的分數(shù)階非線性模型可以更加準確描述橡膠動態(tài)特性。

由于恢復力包含三次方非線性項，系統(tǒng)的動力學響應更加復雜，低頻區(qū)甚至存在超諧波共振。零部件的疲勞破壞與不同頻率的幅值相關(guān)，忽略低頻的幅值會對零部件的疲勞壽命設計帶來誤差[19]。對于隔振系統(tǒng)而言，準確預測每個頻帶區(qū)間的動力學響應是有必要的。精密儀器設備經(jīng)常處于低頻微幅的振動環(huán)境中，金屬橡膠隔振器是一種廣泛應用于微振動隔離的裝置[20-22]。隨著激勵頻率改變，系統(tǒng)的動力學響應也會發(fā)生變化，甚至出現(xiàn)分岔和混沌，進而影響了系統(tǒng)的穩(wěn)定性，對系統(tǒng)的隔振性能以及疲勞壽命預測帶來影響。非線性系統(tǒng)還存在周期運動多樣性，運動轉(zhuǎn)遷過程更加復雜，甚至存在多態(tài)共存現(xiàn)象，這影響了低頻區(qū)隔振系統(tǒng)動力學響應的準確預測。針對整數(shù)階非線性系統(tǒng)的動態(tài)特性研究已相對成熟，而分數(shù)階非線性系統(tǒng)的周期運動多樣性及轉(zhuǎn)遷規(guī)律有待揭示[23-26]。

為了揭示分數(shù)階非線性隔振系統(tǒng)的動力學響應，本文對系統(tǒng)的超諧波共振和周期運動多樣性研究分析。首先對分數(shù)階項進行化簡處理，其次采用高階諧波平衡法求解系統(tǒng)的穩(wěn)態(tài)響應，數(shù)值仿真系統(tǒng)的動力學性能并對近似解析結(jié)果進行比較。接著建立兩個Poincaré映射描述系統(tǒng)周期運動的多樣性及周期運動的轉(zhuǎn)遷規(guī)律，采用Floquet理論對分岔類型加以判定。最后對超諧波共振和周期運動多樣性之間的關(guān)系分析研究，進一步揭示系統(tǒng)參數(shù)對幅頻特性、分岔及混沌的影響。

1 模型建立及響應求解

1.1 分數(shù)階非線性Zener模型建立及分數(shù)階項處理

采用分數(shù)階非線性Zener模型描述橡膠的黏彈性如圖1所示，其中M為系統(tǒng)的質(zhì)量；F為外激勵幅值；X1為質(zhì)量塊的位移；X2為節(jié)點的位移；Ω為外激勵頻率；K1為線性彈性恢復力剛度系數(shù)；C為黏性阻尼系數(shù)；FK為三次非線性彈性恢復力；D為分數(shù)階微分項；K為分數(shù)階項的系數(shù)。引入外激勵幅值F的對照參數(shù)進Fs，進行如下變量代換

圖1 分數(shù)階非線性Zener模型Fig.1 Fractional nonlinear Zener model

分數(shù)階非線性Zener模型無量綱微分方程為

(1)

(2)

Γ(y+1)=yΓ(y)

(3)

式中,Γ(n)為Gamma函數(shù)，滿足式(3)。設質(zhì)量塊的穩(wěn)態(tài)響應為

(4)

將式(4)代入式(2)并引入文獻[24]中的公式

(5a)

(5b)

取分數(shù)階項的一階近似化簡可得

(6)

式(6)所得分數(shù)階項不僅具有阻尼作用也具有剛度作用與很多文獻結(jié)果相同，其中文獻[16]通過參數(shù)識別驗證了分數(shù)階項采用一階三角函數(shù)表示在實際工程應用的合理性。本文主要研究隔振系統(tǒng)高次諧波幅值相比基波幅值是小量的情況，進而分數(shù)階化簡忽略高階項，只取一階近似。

1.2 等效線性剛度系數(shù)和等效阻尼系數(shù)分析

將分數(shù)階項化簡所得一階三角函數(shù)式(6)代入到式(1)中，系統(tǒng)的無量綱微分方程變?yōu)?/p>

βx1(t)]=fcos(ωt)

(7)

系統(tǒng)等效線性剛度系數(shù)Keq和等效阻尼系數(shù)Ceq為

(8)

(9)

由式(8)和式(9)可知，分數(shù)階的系數(shù)λ、分數(shù)階的階數(shù)p對等效線性剛度系數(shù)Keq和等效阻尼系數(shù)Ceq都有一定影響。下面對系統(tǒng)的等效線性剛度系數(shù)Keq和等效阻尼系數(shù)Ceq進行研究。選取系統(tǒng)結(jié)構(gòu)參數(shù)ξ=0.1，μk=2，ε=0.2，f=5，取分數(shù)階的階數(shù)p=0.5。由圖2(a)可以看出分數(shù)階系數(shù)λ越大，等效線性剛度系數(shù)Keq越大，且隨頻率ω的增大而增大。分數(shù)階系數(shù)λ越小，等效阻尼系數(shù)Ceq越小，且隨頻率ω的增大在低頻區(qū)快速減小，離開低頻區(qū)等效阻尼系數(shù)Ceq無明顯變化。其次取分數(shù)階系數(shù)λ=0.01，由圖2(b)可以看出頻率較小時，分數(shù)階項的階數(shù)p越小，等效線性剛度系數(shù)Keq越大；而頻率ω較大時，分數(shù)階的階數(shù)p越小，等效線性剛度系數(shù)Keq越小；等效線性剛度系數(shù)Keq隨著頻率ω增大而增大。在超低頻區(qū)，分數(shù)階的階數(shù)p越小，等效阻尼系數(shù)Ceq越大；離開超低頻區(qū)，分數(shù)階的階數(shù)p越大，等效阻尼系數(shù)Ceq越大。在低頻區(qū)等效阻尼系數(shù)Ceq隨著頻率ω的增大迅速減小，隨著頻率ω繼續(xù)增大離開低頻區(qū)，等效阻尼系數(shù)Ceq無明顯變化。

(a) 系數(shù)λ與等效線性剛度系數(shù)和等效阻尼系數(shù)的關(guān)系

(b) 階數(shù)p與等效線性剛度系數(shù)和等效阻尼系數(shù)的關(guān)系圖2 分數(shù)階項與等效線性剛度系數(shù)和阻尼系數(shù)的關(guān)系Fig.2 Relationship between fractional order term and equivalent linear stiffness and damping

不同分數(shù)階的系數(shù)λ和分數(shù)階的階數(shù)p會改變系統(tǒng)的等效阻尼系數(shù)Ceq和等效線性剛度系數(shù)Keq，進而會改變系統(tǒng)的固有頻率。

1.3 諧波平衡法求解系統(tǒng)穩(wěn)態(tài)響應

設節(jié)點的穩(wěn)態(tài)響應為

(10)

其中(2n-1)表示高階諧波的階次。當n=1時根據(jù)式(4)和式(10)可得質(zhì)量塊和節(jié)點的一階諧波響應為

x1=A1cos(ωt)+B1sin(ωt)=A11sin(ωt+θ1)

(11a)

x2=a1cos(ωt)+b1sin(ωt)=a11sin(ωt+φ1)

(11b)

將式(11)代入到式(1)中進行一階諧波平衡，可得質(zhì)量塊與節(jié)點幅值和相位的關(guān)系

(12)

將式(12)代入到式(11a)中可得到質(zhì)量塊的一階穩(wěn)態(tài)響應

(13)

將式(13)和式(11b)代入式(7)可得

(14)

(15)

進而可以求解出節(jié)點的一階幅值與相位

根據(jù)以上方法獲取的水體信息如圖1所示，同時根據(jù)Landsat OLI成像的準確時間，借助于磨盤山水庫同一時間的水文觀測資料、水庫水位和流量等水文數(shù)據(jù)，基于庫容曲線計算出當時的庫區(qū)水面面積值為22.012 km2，并將不同方法所獲取的水庫面積與實測值進行對比，如表1所示。

(16)

(17)

系統(tǒng)的高次超諧波響應需先求解節(jié)點處的各階諧波幅值系數(shù)，再選取時間歷程圖中節(jié)點的最大位移作為幅值，進而得到質(zhì)量塊的高次響應。當n=2時采用三階諧波平衡法設節(jié)點處的響應為

(18)

對式(18)進行三階諧波平衡可得如下方程組。

(19a)

(19b)

(19c)

36εω3a1b1b3ξ3-30εω2a1a3b1ξ2-3εωa1b1b3ξ+

(19d)

系統(tǒng)參數(shù)確定的條件下借助計算機求解并對結(jié)果進行數(shù)據(jù)處理，即可得到節(jié)點幅值系數(shù)a1、b1、a3、b3的值[17]，再將節(jié)點幅值系數(shù)代入式(1)可得質(zhì)量塊的幅值系數(shù)

(20)

進而得到質(zhì)量塊的三次穩(wěn)態(tài)響應。

(21)

根據(jù)三階諧波平衡法的求解過程以此類推可得系統(tǒng)的高次諧波響應。

1.4 數(shù)值方法對比

隨機選取一組系統(tǒng)參數(shù)：μk=0.5，f=5，ξ=0.1，ε=0.1，λ=0.1。采用四階Runge-Kutta法對系統(tǒng)數(shù)值求解，通過動力學仿真軟件UM對隔振系統(tǒng)進行虛擬實驗仿真，不同方法求解系統(tǒng)動力學響應結(jié)果接近，如圖3所示。

(a) 質(zhì)量塊瞬態(tài)響應

(b) 質(zhì)量塊穩(wěn)態(tài)響應圖3 質(zhì)量塊時間歷程圖Fig.3 Time history diagram of mass block

一階諧波平衡法只能求解系統(tǒng)的主共振，無法求解系統(tǒng)的超諧波共振，為進一步研究系統(tǒng)在低頻區(qū)的動力學響應，采用三階諧波平衡法對系統(tǒng)進行求解進而得到超諧波共振如圖4所示。

圖4 不同方法求解系統(tǒng)的幅頻響應曲線Fig.4 Amplitude-frequency response curve by different methods

2 系統(tǒng)運動狀態(tài)與分岔類型

3 超諧波共振下的周期運動多樣性

為方便描述系統(tǒng)運動規(guī)律，定義符號T-N來描述系統(tǒng)的運動狀態(tài)，其中T為系統(tǒng)運動的周期數(shù)，N為周期運動內(nèi)簡諧振動的次數(shù)。為避免混沌出現(xiàn)，選取阻尼系數(shù)相對較大的基準參數(shù)：μk=0.15，ε=0.6，p=0.5，λ=0.1，f=5，ξ=0.5。選取n=5采用高階諧波平衡法可得系統(tǒng)質(zhì)量塊的幅頻曲線如圖5所示。

(a) 主共振及超諧波共振

(b) 超諧波共振局部放大圖圖5 幅頻特性曲線圖Fig.5 Amplitude-frequency response curve

從圖5中可見高次超諧波共振幅值相比較主共振幅值非常小，但是在掃頻的過程中高次超諧波也存在跳躍現(xiàn)象如圖5(b)所示，系統(tǒng)在低頻區(qū)也存在多態(tài)共存現(xiàn)象且動力學行為更加復雜。隨著激勵頻率改變，系統(tǒng)的動力學響應出現(xiàn)分岔和混沌，這影響了隔振系統(tǒng)的穩(wěn)定性，對系統(tǒng)的動力學響應預測帶來影響。隨著頻率ω的降低，數(shù)值模擬系統(tǒng)的運動狀態(tài)規(guī)律：1-1→1-2→1-1→1-3→1-5→1-7→1-9，如圖6所示。

(a) 定相位面

(b) 定極大位移面圖6 系統(tǒng)質(zhì)量塊的運動狀態(tài)Fig.6 Motion state of system mass block

對比超諧波共振的階次可見超諧波共振的階次越高，周期運動內(nèi)的簡諧振動次數(shù)越多，周期運動內(nèi)簡諧振動的次數(shù)與超諧波共振的階次近似一致。隨著超諧波次數(shù)的增加，幅頻曲線跳躍的距離減小且跳躍方向與分岔圖中SNB突變的方向一致。由于設解的形式均為奇數(shù)次諧波，無論諧波平衡法設其解的階次為多大，該方法都具有一定局限性。從圖7相圖可以看到在ω<1的低頻區(qū)，隨著頻率ω減小系統(tǒng)始終保持周期一運動，但是周期內(nèi)的簡諧振動次數(shù)增加。

(a) ω=1的Poincaré截面圖和相圖(1-1運動)

(b) ω=0.696 1的Poincaré截面圖和相圖(1-3運動)

(c) ω=0.440 6的Poincaré截面圖和相圖(1-5運動)

(d) ω=0.328 6的Poincaré截面圖和相圖(1-7運動)

(e) ω=0.230 6的Poincaré截面圖和相圖(1-9運動)圖7 系統(tǒng)的Poincaré截面圖和相圖Fig.7 Poincaré section diagram and phase diagram of the system

由于Duffing系統(tǒng)的恢復力為奇函數(shù)，這也導致系統(tǒng)的運動狀態(tài)在相圖中也呈現(xiàn)為反對稱的運動軌跡，運動狀態(tài)及運動轉(zhuǎn)遷過程更加復雜。

4 分數(shù)階項對超諧共振及周期運動轉(zhuǎn)遷的影響

下面主要研究分數(shù)階微分項的系數(shù)λ和階數(shù)p對系統(tǒng)超諧波共振和周期運動轉(zhuǎn)遷的影響，選取n=5采用高階諧波平衡法求解系統(tǒng)響應，數(shù)值仿真系統(tǒng)動力學響應。

首先分析分數(shù)階系數(shù)λ對超諧波共振的影響，選取基準參數(shù)為：μk=0.15，ε=0.6，f=5，ξ=0.01，取分數(shù)階的階數(shù)p=0.5，分數(shù)階系數(shù)分別取λ=0.1、λ=0.2、λ=0.3，系統(tǒng)在低頻區(qū)的超諧波共振如圖8(a)所示。由于分數(shù)階項不僅具有阻尼作用也具有剛度作用，進而隨著分數(shù)階系數(shù)λ增大，超諧波共振的幅值與彎曲程度減小。當分數(shù)階系數(shù)繼續(xù)增大到λ=0.3時，9次超諧波共振消失，可見較大的分數(shù)階系數(shù)λ不僅可以有效降低超諧波共振的幅值，還可以抑制高次超諧波共振的出現(xiàn)。

(a) 不同分數(shù)階的系數(shù)λ

(b) 不同分數(shù)階的階數(shù)p圖8 超諧波共振幅頻曲線圖Fig.8 Super harmonic resonance amplitude frequency characteristic curve

其次分析分數(shù)階的不同階數(shù)p對系統(tǒng)超諧波共振的影響。選取基本參數(shù)為μk=0.15，ε=0.6，f=5，ξ=0.01，分數(shù)階系數(shù)取λ=0.1，分數(shù)階的階數(shù)分別取p=0.2、p=0.5、p=0.8，系統(tǒng)在低頻區(qū)的超諧波共振如圖8(b)所示。從圖8(b)可見在頻率ω<0.3的局部放大圖中，分數(shù)階的階數(shù)p越小，7、9次超諧波共振的峰值越大；當頻率ω>0.3時，3、5次超諧波共振的峰值隨階數(shù)p的增加而增高。

接著研究分數(shù)階系數(shù)λ對系統(tǒng)周期運動轉(zhuǎn)遷的影響。數(shù)值仿真得到不同分數(shù)階系數(shù)λ所對應系統(tǒng)質(zhì)量塊的分岔圖如圖9所示，虛線表示系統(tǒng)發(fā)生鞍結(jié)分岔SNB。從圖9可見分數(shù)階系數(shù)λ由0.1增加到0.2時，分岔圖中的SNB的個數(shù)由4減小為1并且3次超諧波和5次超諧波轉(zhuǎn)遷過程中的混沌消失；隨著分數(shù)階系數(shù)λ繼續(xù)增加為0.3時，鞍結(jié)分岔SNB也消失。可見分數(shù)階系數(shù)λ的增大還可有效避免系統(tǒng)出現(xiàn)分岔和混沌。

(a) λ=0.1的頻率正掃分岔圖和頻率反掃分岔圖

(b) λ=0.2的頻率正掃分岔圖和頻率反掃分岔圖

(c) λ=0.3的頻率正掃分岔圖和頻率反掃分岔圖圖9 不同分數(shù)階系數(shù)對應質(zhì)量塊的分岔圖Fig.9 Bifurcation diagram of mass block corresponding to different fractional order coefficients

其次研究分數(shù)階的不同階數(shù)p對系統(tǒng)周期運動轉(zhuǎn)遷規(guī)律的影響如圖10所示，虛線表示系統(tǒng)發(fā)生鞍結(jié)分岔SNB。當分數(shù)階的階數(shù)p由0.2增加為0.5時，SNB的個數(shù)由5減小為4，在頻率ω=0.234 2的SNB消失，在頻率ω=0.31處的混沌消失；隨著分數(shù)階的階數(shù)p繼續(xù)增加為0.8時，SNB的個數(shù)由4又增加為5，在頻率ω=0.234 2處再次出現(xiàn)SNB，可見在本節(jié)基準參數(shù)下為有效抑制在超低頻區(qū)出現(xiàn)分岔和混沌，分數(shù)階的階數(shù)p應取值適中。

(a) p=0.2的頻率正掃分岔圖和頻率反掃分岔圖

(b) p=0.5的頻率正掃分岔圖和頻率反掃分岔圖

對比幅頻曲線和分岔圖可見在相鄰次數(shù)超諧波共振的轉(zhuǎn)遷過渡區(qū)存在混沌，并且存在多態(tài)共存現(xiàn)象。為進一步確定多態(tài)共存下的運動狀態(tài)及轉(zhuǎn)遷規(guī)律，在本節(jié)選取基準參數(shù)下，選取分數(shù)階的系數(shù)λ=0.1和分數(shù)階的階數(shù)p=0.5，系統(tǒng)多態(tài)共存區(qū)域主要分布在主共振和3次超諧波共振轉(zhuǎn)遷過渡區(qū)，以及3次超諧波共振和5次超諧波共振轉(zhuǎn)遷過渡區(qū)，系統(tǒng)運動狀態(tài)的多樣性如圖11所示。

圖11 多態(tài)共存及其相鄰區(qū)域分岔圖Fig.11 Multi state coexistence and its adjacent region bifurcation diagram

由于系統(tǒng)存在SNB、倍周期分岔(PDB)、叉式分岔(PFB)以及邊界激變(BC)等多種分岔，在這些分岔的誘導下甚至會出現(xiàn)混沌，并且該過程不可逆，具體轉(zhuǎn)遷規(guī)律如圖12所示。

圖12 多態(tài)共存及相鄰區(qū)域運動轉(zhuǎn)遷規(guī)律Fig.12 Multi state coexistence and movement transition law of adjacent regions

在橡膠工業(yè)中，可根據(jù)上述規(guī)律選擇恰當?shù)南到y(tǒng)參數(shù)，進而避免隔振系統(tǒng)動力學響應出現(xiàn)非線性跳躍現(xiàn)象以及分岔和混沌，使橡膠隔振系統(tǒng)的動力學性能更佳。

5 結(jié) 論

本文采用分數(shù)階非線性Zener模型描述黏彈性隔振系統(tǒng)，采用高階諧波平衡法求解了系統(tǒng)的穩(wěn)態(tài)響應并結(jié)合多種方法對結(jié)果進行比較，闡述了低頻區(qū)系統(tǒng)幅頻特性與周期運動多樣性之間的關(guān)系，研究了分數(shù)階項對中低頻范圍內(nèi)隔振系統(tǒng)的分岔、混沌和多態(tài)共存等復雜非線性動力學行為的影響，得出以下結(jié)論：

(1) 針對系統(tǒng)在低頻區(qū)域的超諧波共振，分岔圖中出現(xiàn)鞍結(jié)分岔SNB對應的頻率與幅頻曲線發(fā)生跳躍的頻率一致，鞍結(jié)分岔SNB突變的距離隨諧波次數(shù)增加而減小。

(2) 分數(shù)階系數(shù)λ的增大不僅可有效降低超諧波共振的峰值，還可有效避免系統(tǒng)發(fā)生分岔和混沌。

(3) 相鄰次數(shù)超諧波共振轉(zhuǎn)遷過程中存在多態(tài)共存的現(xiàn)象，轉(zhuǎn)遷過程中相繼出現(xiàn)PFB、SNB、PDF和BC等多種分岔類型，在分岔誘導下甚至會出現(xiàn)混沌且該過程不可逆。

上述研究結(jié)果與方法可為高次諧波為小量的隔振系統(tǒng)在微振動下的動力學響應預測提供一定理論依據(jù)。