分數階非線性隔振系統的超諧波共振與周期運動轉遷規律分析

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

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

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

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

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

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

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

1 模型建立及響應求解

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

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

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

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

(1)

(2)

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

(3)

式中,Γ(n)為Gamma函數，滿足式(3)。設質量塊的穩態響應為

(4)

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

(5a)

(5b)

取分數階項的一階近似化簡可得

(6)

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

1.2 等效線性剛度系數和等效阻尼系數分析

將分數階項化簡所得一階三角函數式(6)代入到式(1)中，系統的無量綱微分方程變為

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

(7)

系統等效線性剛度系數Keq和等效阻尼系數Ceq為

(8)

(9)

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

(a) 系數λ與等效線性剛度系數和等效阻尼系數的關系

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

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

1.3 諧波平衡法求解系統穩態響應

設節點的穩態響應為

(10)

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

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

(11a)

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

(11b)

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

(12)

將式(12)代入到式(11a)中可得到質量塊的一階穩態響應

(13)

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

(14)

(15)

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

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

(16)

(17)

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

(18)

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

(19a)

(19b)

(19c)

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

(19d)

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

(20)

進而得到質量塊的三次穩態響應。

(21)

根據三階諧波平衡法的求解過程以此類推可得系統的高次諧波響應。

1.4 數值方法對比

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

(a) 質量塊瞬態響應

(b) 質量塊穩態響應圖3 質量塊時間歷程圖Fig.3 Time history diagram of mass block

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

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

2 系統運動狀態與分岔類型

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

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

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

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

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

(a) 定相位面

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

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

(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 系統的Poincaré截面圖和相圖Fig.7 Poincaré section diagram and phase diagram of the system

由于Duffing系統的恢復力為奇函數，這也導致系統的運動狀態在相圖中也呈現為反對稱的運動軌跡，運動狀態及運動轉遷過程更加復雜。

4 分數階項對超諧共振及周期運動轉遷的影響

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

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

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

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

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

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

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

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

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

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

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

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

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

圖11 多態共存及其相鄰區域分岔圖Fig.11 Multi state coexistence and its adjacent region bifurcation diagram

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

圖12 多態共存及相鄰區域運動轉遷規律Fig.12 Multi state coexistence and movement transition law of adjacent regions

在橡膠工業中，可根據上述規律選擇恰當的系統參數，進而避免隔振系統動力學響應出現非線性跳躍現象以及分岔和混沌，使橡膠隔振系統的動力學性能更佳。

5 結 論

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

(1) 針對系統在低頻區域的超諧波共振，分岔圖中出現鞍結分岔SNB對應的頻率與幅頻曲線發生跳躍的頻率一致，鞍結分岔SNB突變的距離隨諧波次數增加而減小。

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

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

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