雙頻激勵下含分數階非線性汽車懸架系統的混沌研究

常宇健，孫亞婷，陳恩利，李韶華，邢武策

（1.石家莊鐵道大學電氣與電子工程學院，河北石家莊050043；2.石家莊鐵道大學省部共建交通工程結構力學行為與系統安全國家重點實驗室，河北石家莊050043）

引言

分數階微積分作為重要的數學分支，于1695年德國科學家Leibniz和法國數學家L′Hopital在探討1/2階導數時首次被提出［1］。然而，由于缺乏應用背景支撐等多方面原因，它長期以來并沒有得到較多的關注和研究。隨著20世紀70年代以來對分形和各種復雜系統的深入研究，分數階微積分理論及其應用開始受到廣泛關注，很多學者對分數階微積分的基本特性進行研究，在基礎理論方面取得了很大進展［2-6］。

進入21世紀以來，分數階微積分建模方法和理論在復雜黏彈性材料力學本構關系、反常擴散、高能物理等諸多領域有了若干非常成功的應用［7-10］，凸顯了其獨特優勢和不可替代性，所以研究含分數階微積分方程中的典型力學特性和分數階參數對動力系統的影響很有意義，大量學者進行著這方面的研究［11-13］。

車輛懸架減振裝置不僅具有遲滯非線性特性，而且多數阻尼器都具有類黏彈性本構關系，這些黏彈性材料介于彈性和阻尼特性之間，普通的整數階理論無法準確地描述這種材料的本構關系。為此，很多研究學者開始對黏彈性材料采用分數階理論進行描述［14-16］。將分數階理論和非線性理論一起運用到汽車懸架系統中，不僅具有一定的前瞻性，而且由于在懸架系統中加入了分數階微分理論，可以更準確地描述懸架系統的數學模型。因此本文采用含分數階微分項的非線性方程來描述汽車懸架系統。

在外界激勵的作用下，具有遲滯非線性的汽車懸架系統會產生復雜的非線性動力學行為，如分叉、混沌等。劉劍等［17］建立了兩自由度磁流變懸架動力學系統，根據非線性穩定性理論發現了單頻諧波激勵下系統發生混沌的可能性。黃苗玉等［18］針對兩自由度磁流變懸架系統，研究了簡諧路面作用下系統隨激勵頻率幅值變化的分叉特性，并利用相平面圖、龐加萊截面圖等詳細描述了通向混沌振動的路徑。樓京俊等［19］利用Melnikov方法研究多頻激勵下的軟彈簧Duffing系統的混沌動力學，證明激勵頻率數目的增加使系統更容易進入混沌狀態。畢勤勝等［20］討論了參外聯合激勵復合非線性振子的動力學行為，并對其定常解進行了局部分叉分析。楊智勇等［21］建立了雙頻擬周期動態路面激勵函數，并構建了四自由度非線性車輛懸架模型，通過分析系統的龐加萊圖、相位圖等得到了系統發生混沌時的激勵振幅和振動特性。李韶華等［22］研究了具有滯后非線性的汽車懸架在路面擬周期激勵作用下發生受迫振動的混沌運動，揭示出此系統中存在從擬周期運動通向混沌運動的可能性。

目前對汽車懸架的研究主要集中在系統的響應分析及懸架的控制技術，對雙頻激勵非線性系統的動力學行為研究，尤其是針對含分數階的非線性系統在雙頻激勵下的動力學行為研究較少。針對雙頻激勵下的含分數階非線性汽車懸架系統，本文通過計算Melnikov函數得出了系統的混沌閾值，并分析了參數對混沌區域的影響。

1 雙頻激勵下Smale馬蹄變換意義下混沌的解析必要條件

1.1 含分數階非線性懸架系統建模

汽車懸架非線性模型如圖1所示。圖中：z0為路面激勵位移，z1為系統垂直位移。其運動微分方程為

式中m，k1，k3，c1，h和p分別為懸架系統的質量、線性剛度系數、非線性剛度系數、阻尼系數、分數階項的系數和階次。

設：相對位移z=z1-z0，則z1=z+z0，式（1）可以轉化為

由式（2）可知，含分數階非線性懸架模型是一種分數階杜芬系統。擬周期激勵表示為：z0=其中，F1，F2為激勵幅值；ω1，ω2為激勵頻率，且ω1，ω2不可有理通約。令，式（2）可以轉化為

對式（3）進行變量代換，令

則有

1.2 分數階杜芬系統的異宿軌道

將式（3）寫成狀態方程的形式

通過M點和N點的異宿軌道滿足

對式（7）進行整理可得

對式（8）進行分離變量與積分分解，從而得到異宿軌道為

1.3 Melnikov函數確定系統混沌邊界條件

系統的Melnikov函數如下

根據Melnikov理論可得發生混沌的邊界條件為

其中，

對于式（11）中的I3，由于被積函數中含有分數階微分項，不能直接進行積分運算，只能通過分數階微分的定義來進行計算。分數階導數的定義不像常規整數階導數有唯一的定義式，根據不同的研究背景，分數階導數在其發展過程中被給予了多種形式的定義式。目前使用較多的有Riemann-Liouville定義、Caputo定 義 及Grunwald-Letnilov定 義。Riemann-Liouville定義具有良好的數學性質，但在工程中應用受到諸多限制。例如，初值問題在工程中具有重要意義，但常數的分數階導數在Riemann-Liouville定義下卻不為零。Caputo定義在工程應用中，其物理意義更加明確，對于初始條件非零的微分方程，Caputo分數階微分可以降低其計算難度，但Caputo定義下對分數階微積分進行離散計算比較困難。Riemann-Liouville形式的分數階導數方便進行數值運算，但易于離散。本文使用Caputo形式的分數階導數定義法對分數階微分方程進行數值計算。Caputo與Riemann-Liouville分數階微分的定義分別為：

式中Γ(x)為Gamma函數，且Γ(x+1)=xΓ(x)。

通過比較式（12）與（13）分數階微分的兩種定義可知，當x(t)的初值非零時，兩者之間的關系為

計算高精度Caputo分數階微分導數的具體過程為：首先，通過改進的直接遞推算法［23］來對Riemann-Liouville分數階微分進行高精度計算；然后，計算出補償函數f(t)；最后，將補償函數與Riemann-Liouville分數階微分相結合，得到高精度的Caputo分數階微分，計算結果可靠。這樣可以計算出含有分數階微分項的積分，得到I3的值。

將計算得出的I1，I2，I3，I4帶入邊界條件式（11）得到系統發生具有Smale馬蹄意義下混沌的必要條件表達式為

2 含分數階非線性懸架系統參數對混沌閾值的影響

根據式（15）所求出的激勵幅值F與激勵頻率ω1，ω2之 間 的 關 系，取m=240 kg，c1=200 N?s/m，k1=15000 N/m，k3=30000 N/m3，h=1000，p=0.5，可以確定系統混沌的邊界條件，所得混沌邊界曲線如圖2所示。

圖2 混沌邊界曲面Fig.2 Chaotic boundary surface

當F的取值位于混沌邊界曲面之上時，則對于充分小的0＜ε≤1，系統可能發生Smale馬蹄意義下的混沌；若F位于混沌邊界曲面下方，系統做穩定的擬周期運動。分別改變c1，k1，k3，h和p的值，可研究系統各參數對混沌邊界曲面的影響。其中圖3為阻尼系數c1的影響（沿箭頭方向c1取值依次增大，分別為160，180，200，220，240）；圖4為線性剛度系數k1的影響（沿箭頭方向k1取值依次增大，分別為11000，13000，15000，17000，19000）；圖5為非線性剛度系數k3的影響（沿箭頭方向k3取值依次增大，分別為26000，28000，30000，32000，34000）；圖6為分數階項系數h的影響（沿箭頭方向h取值依次增大，分別為600，800，1000，1200，1400）；圖7為分數階項階數的影響（沿箭頭方向p取值分別為0.3，0.1，0，0.5，0.7，0.9，1）。

圖4 k1對混沌邊界曲面的影響Fig.4 Influence of k1 on the chaotic boundary surface

圖5 k3對混沌邊界曲面的影響Fig.5 Influence of k3 on the chaotic boundary surface

圖6 h對混沌邊界曲面的影響Fig.6 Influence of h on the chaotic boundary surface

圖7 p對混沌邊界曲面的影響Fig.7 Influence of p on the chaotic boundary surface

從圖中可以看出：

（1）c1越大，混沌閾值越小，發生混沌的可能性越大，閾值隨c1變化比較均勻。

（2）k1越大，混沌閾值越小，發生混沌的可能性越大，閾值隨k1變化也比較均勻。

（3）k3越大，混沌閾值越小，發生混沌的可能性越大，但非線性項k3的變化對閾值的影響并不十分明顯。

（4）h越大，混沌閾值越大，發生混沌的可能性越小，h的改變對閾值變化有明顯的影響。

（5）p從0變化到1時，混沌閾值先減小后增大，發生混沌的可能性先增大后減小，且階數變化對閾值影響較大。

由圖可知，懸架系統的各個參數對系統的混沌邊界曲線均有影響，其中分數階微分項的階數和系數對懸架系統的混沌邊界曲線影響很大。由此可知，將分數階微分項引入汽車懸架系統，比單純使用整數階對懸架系統進行描述更能準確描述系統的混沌閾值，這對實際工程中懸架參數的設計和選擇具有一定的參考意義和理論價值。

3 數值仿真驗證及分析

為了驗證邊界的正確性，本文隨機研究了擬周期頻率ω1=7.3 rad/s，ω2=20 rad/s及ω1=8.1 rad/s，ω2=21 rad/s時的情況。首先分數階項使用冪級數展開法進行數值運算，根據Grunwald-Letnilov形式分數階導數定義對其進行近似，并將代數方程離散化即可。然后根據系統在雙頻激勵下的時間歷程圖、頻譜圖、相圖及龐加萊截面圖進行判斷，最后根據最大Lyapunov指數進一步驗證。

3.1 不同激勵頻率及幅值下的數值仿真運動狀態結果研究

進行數值仿真的兩組數據在混沌邊界曲面的實際位置如圖8所示。

圖8 各點在混沌邊界曲面的實際位置Fig.8 Position of each point on the chaotic boundary surface

通過數值仿真分析可以得出系統在混沌曲面上下各點的時間歷程圖、頻譜圖、相圖及龐加萊截面圖如圖9-17所示。表1為各點取值、位置及運動狀態。根據數值仿真結果判斷系統的運動形式，進而驗證混沌邊界曲面的正確性。下面分別對周期運動區域及可能發生混沌的區域選取兩組不同激勵頻率下的點進行比較驗證。

表1 各點位置及運動狀態Tab.1 Position and motion state of each point

從圖9-13可以得出：圖9和10其時間歷程圖為擬周期運動，頻譜圖上只有兩個與外界激勵頻率相等的主頻率，相圖表現為穩定的相軌線，龐加萊截面圖類似一個圓環。因此，可以判斷系統在此外界激勵下做穩定的擬周期運動。A點位于混沌閾值曲面的下方，為擬周期運動。B點雖然位于混沌閾值曲面上方，但解析求得的閾值曲面上方為發生混沌的必要條件，并非充要條件，因此B點仍為擬周期運動。

圖9 A點擬周期運動(ω1=7.3 rad/s,ω2=20 rad/s,F=0.02 m)Fig.9 Quasi periodic motion of point A(ω1=7.3 rad/s,ω2=20 rad/s,F=0.02 m)

圖10 B點 擬 周 期 運 動(ω1=7.3 rad/s,ω2=20 rad/s,F=0.05 m)Fig.10 Quasi periodic motion of point B(ω1=7.3 rad/s,ω2=20 rad/s,F=0.05 m)

C點位于混沌閾值曲面上方，由圖11可以看出時間歷程圖基本沒有規律，頻譜圖上出現多個離散的譜線，相圖開始趨于混亂，龐加萊截面圖的圓環出現破裂。因此，初步判斷系統在此激勵下系統開始出現混沌運動，但比較接近于擬周期運動。

圖11 C點混沌運動(ω1=7.3 rad/s，ω2=20 rad/s，F=0.1 m)Fig.11 Chaotic motion of point C(ω1=7.3 rad/s，ω2=20 rad/s，F=0.1 m)

隨著外界激勵幅值逐漸增大，從圖12和13可以看出，時間歷程圖越來越混亂無序，頻譜圖出現多個離散譜線，且在低頻部分表現為連續混亂的頻譜，相圖軌線越來越混亂，龐加萊截面圖圓環完全破裂，呈現為毫無規律的一些離散點，可以判斷系統呈現出完全的混沌運動。

圖12 D點混沌運動(ω1=7.3 rad/s，ω2=20 rad/s，F=0.12 m)Fig.12 Chaotic motion of point D(ω1=7.3 rad/s，ω2=20 rad/s，F=0.12 m)

從由圖14-17可以得出：圖14和15的時間歷程圖為規則的擬周期運動，頻譜圖只有兩條與外界激勵頻率相等的主頻率，相圖均表現為穩定的相軌線，龐加萊截面圖為一個整齊的圓環。因此，可以判斷系統在此外界激勵下做穩定的擬周期運動。

圖13 E點 混 沌 運 動（ω1=7.3 rad/s，ω2=20 rad/s，F=0.13 m）Fig.13 Chaotic motion of point E（ω1=7.3 rad/s，ω2=20 rad/s，F=0.13 m）

圖14 F點 擬 周 期 運 動（ω1=8.1 rad/s,ω2=21 rad/s,F=0.02 m）Fig.14 Quasi periodic motion of point F（ω1=8.1 rad/s,ω2=21 rad/s,F=0.02 m）

圖15 G點 擬 周 期 運 動（ω1=8.1 rad/s,ω2=21 rad/s,F=0.03 m）Fig.15 Quasi periodic motion of point G（ω1=8.1 rad/s,ω2=21 rad/s,F=0.03 m）

從圖16可以看出，系統的時間歷程圖表現為較為整齊規則的擬周期運動，頻譜圖出現多條離散譜，相圖開始趨于混亂無序，龐加萊截面圖中圓環開始出現破裂。初步判斷系統在此外界激勵下做擬周期運動，但已經比較接近混沌運動。

圖16 H點 擬 周 期 運 動（ω1=8.1 rad/s,ω2=21 rad/s,F=0.13 m）Fig.16 Quasi periodic motion of point H（ω1=8.1 rad/s,ω2=21 rad/s,F=0.13 m）

從圖17可以看出，系統的時間歷程圖呈混亂無序的混沌運動，頻譜圖低頻段出現連續頻譜，相圖出現混亂軌線，龐加萊截面圖為一些雜亂的離散點。可以判斷在此外界激勵下，系統做混沌運動。

圖17 I點混沌運動（ω1=8.1 rad/s,ω2=21 rad/s,F=0.15 m）Fig.17 Chaotic motion of point I（ω1=8.1 rad/s,ω2=21 rad/s,F=0.15 m）

3.2 最大Lyapunov指數原理及結果

含分數階非線性系統的振動狀態具有非常復雜的動力學特性，可以利用最大Lyapunov指數表征系統在相空間中相鄰軌道的發散率，若最大Lyapunov指數為正值，則表示兩條相鄰軌線隨時間變化成指數率增長，系統運動具有混沌狀態；若最大Lyapunov指數為負，則表示系統對初值不敏感，系統狀態收斂到平衡點，此時系統具有擬周期運動特性。為驗證上述混沌邊界曲面的正確性，使用MATLAB軟件對含分數階非線性汽車懸架系統進行數值仿真分析。

進行數值仿真分析時，對最大Lyapunov指數使用Wolf重構算法，仿真程序中選取的初始條件為初始位移為-0.5 m，初始速度為0.5858 m/s，分數階項初始值為0.1，步長為計算A-I各點最大Lyapunov指數及運動狀態如表2所示。所有出現混沌的點都位于邊界曲面上方，據此可以進一步判斷其運動狀態的正確性。根據時間歷程圖、頻譜圖、相圖及龐加萊截面圖判斷所得各點運動狀態與根據最大Lyapunov指數判斷結果一致，因此可以驗證所得混沌邊界曲面正確。

表2 各點最大Lyapunov指數及運動狀態Tab.2 Maximum Lyapunov exponent and motion state of each point

4 結論

本文針對含分數階非線性特性的1/4汽車懸架系統，研究雙頻激勵下系統的動力學響應。

（1）通過構建含分數階1/4汽車非線性懸架系統，應用Melnikov解析方法得到系統發生混沌的解析必要條件。

（2）采用數值仿真得到汽車系統響應特征，并結合最大Lyapunov指數分析，驗證了含分數階非線性懸架系統在一定條件下存在混沌運動，且隨著激勵幅值的增加，懸架系統運動狀態為擬周期運動到混沌運動。通過系統最大Lyapunov指數和系統響應驗證了混沌閾值曲面的正確性。

（3）分析表明系統的剛度、阻尼及分數階項系數和階次均對混沌邊界曲面有一定影響，且分數階項參數對系統的影響較大，分析結果可以為懸架系統參數的選擇及混沌行為的研究提供借鑒意義。