基于分?jǐn)?shù)階微積分的油氣懸架建模與試驗(yàn)分析

孫會來，金 純，張文明，李 昊，田海勇

（1．北京科技大學(xué) 機(jī)械工程學(xué)院，北京 100083；2．中國北車股份有限公司，北京 100083）

油氣懸架具有良好的非線性阻尼特性，結(jié)構(gòu)簡單，性能可靠，其特性對車輛行駛平順性和操縱穩(wěn)定性有重要影響，在大客車、載重貨車以及工程機(jī)械車輛上被廣泛采用，特別是在載重自卸車領(lǐng)域應(yīng)用更加普遍。針對油氣懸架建模及其運(yùn)動特性國內(nèi)外學(xué)者做出了大量研究，Concordia大學(xué)的研究人員［1］先后建立了單筒單氣室油氣懸架缸的參數(shù)化計(jì)算模型和非線性模型，采用基于能量的頻率相關(guān)等效線性化技術(shù)，其等效線性參數(shù)化模型被廣泛采用；Gao等［2］對重型車輛主動油氣懸架進(jìn)行了建模及仿真；孫濤等［3］建立了油氣懸架長通孔紊流阻尼模型；楊波等［4］基于懸架非線性阻尼對雙氣室懸架進(jìn)行了建模研究。在上述研究中建立的油氣懸架數(shù)學(xué)模型在考慮油液可壓縮的特性時，多認(rèn)為油液含氣以后體積彈性模量發(fā)生變化。且同時做了氣體的質(zhì)量保持不變，無泄漏，無形變等多項(xiàng)假設(shè)［5－8］。而黃夏旭等［9］研究又表明氣體溶解效應(yīng)和油液可壓縮性在油氣懸架設(shè)計(jì)研究過程中不可忽略。

懸架系統(tǒng)的實(shí)際工作要更為復(fù)雜，不僅與本身參數(shù)有關(guān)，同時受到各處零部件相對運(yùn)動位移、溫度、磨損情況、油液黏度及其在油氣懸架內(nèi)流動及可壓縮特性的影響。以上諸多研究中總是側(cè)重于某幾項(xiàng)參數(shù)的變化及對懸架系統(tǒng)特性的影響，或者是通過增加計(jì)算項(xiàng)來保證計(jì)算精度，同時假設(shè)過多計(jì)算復(fù)雜難以從總體實(shí)現(xiàn)準(zhǔn)確快速描述。

傳統(tǒng)整數(shù)階導(dǎo)數(shù)模型需要通過引入多個導(dǎo)數(shù)項(xiàng)和材料參數(shù)，本文提出根據(jù)油氣懸架多相介質(zhì)的力學(xué)特點(diǎn)，引入分?jǐn)?shù)階微積分理論對運(yùn)動微分方程中的阻尼損耗項(xiàng)進(jìn)行修正建模并分析其等效黏彈性，建立某型工程車輛油氣懸架的振動分?jǐn)?shù)階振動方程并進(jìn)行數(shù)值求解，探討在該模型及不同工況下懸架的振動特性。

1 油氣懸架結(jié)構(gòu)及原理

單氣室油氣懸架主要由缸筒以及活塞桿和活塞組件組成，油氣懸架上端設(shè)有氣室或儲能器。油缸氣室與缸筒上端通過浮動活塞隔開，整個懸架缸內(nèi)形成工作腔和環(huán)形腔，活塞桿壁上設(shè)有阻尼孔和單向閥。工作腔上部充入惰性氣體，下部及環(huán)形腔充入油液。圖1是根據(jù)上述描繪簡化后的油氣懸架物理模型。

圖1 單氣室油氣懸架原理Fig．1 Single-chamber pneumatic suspension principle

圖2 油氣懸架實(shí)車測試Fig．2 Pneumatic suspension substantial vehicle test

當(dāng)懸架處在壓縮行程時，活塞桿向上運(yùn)動，工作腔惰性氣體體積減小，壓力升高，油液受到壓縮經(jīng)單向閥和阻尼孔流入環(huán)形腔，使氣室容積減小，氣壓力升高。在這此過程單向閥及阻尼孔會同時連通工作腔與環(huán)形腔，油液流量大，產(chǎn)生的油液阻尼力比較小，主要由氣室內(nèi)氣體受到壓縮產(chǎn)生彈性作用來抑制活塞桿的向上運(yùn)動。當(dāng)懸架處在伸張行程時，活塞桿向下運(yùn)動，工作腔的體積增大進(jìn)而壓力減小，環(huán)形腔中的油液流向工作腔，此時單向閥關(guān)閉，油液只流經(jīng)阻尼孔，流速較高，產(chǎn)生較大的阻尼力，抑制了活塞桿的向下運(yùn)動，產(chǎn)生的阻尼力較大，有較強(qiáng)衰減振動的作用。工程車輛油氣懸架安裝如圖2所示，油氣彈簧懸架與輪胎組成振動系統(tǒng)起衰減振動、減緩沖擊的作用。

經(jīng)過如上分析及在試驗(yàn)中發(fā)現(xiàn)油液黏度及油液在油氣懸架內(nèi)流動、氣體在油液中的溶解率的隨溫度和壓力而變化等這些都不同程度的影響到了油氣懸架特性。油氣懸架阻尼特性等同于純黏性體的計(jì)算精度有待提高，往往需要涉及過多參數(shù)。油氣懸架特性和外力大小，溫度的改變、力的作用時間及加載歷史都有關(guān)系，這些符合黏彈性物質(zhì)的力學(xué)特點(diǎn)。黏彈性材料的力學(xué)性質(zhì)既表現(xiàn)出彈性，同時又有黏性特性，介于理想彈性體和牛頓流體之間，其應(yīng)力－應(yīng)變響應(yīng)依賴于時間和應(yīng)變率有記憶性，與荷載和變形歷史有關(guān)。而分?jǐn)?shù)階黏彈性微積分模型僅用較少的參數(shù)便能準(zhǔn)確地描述大量復(fù)雜黏彈性材料廣泛頻率范圍內(nèi)的動力學(xué)特性［10］，能較準(zhǔn)確地描述黏彈性材料的動力學(xué)特性。

基于以上思想本文將油氣懸架等效為黏彈性系統(tǒng)進(jìn)行建模，引入分?jǐn)?shù)階模型描述阻尼損耗項(xiàng)對油氣懸架系統(tǒng)進(jìn)行分析，并進(jìn)一步確定分?jǐn)?shù)階次，觀察在不同激勵下分?jǐn)?shù)階建模與整數(shù)階區(qū)別。

2 油氣懸架分?jǐn)?shù)階模型

2．1 分?jǐn)?shù)階微積分定義

分?jǐn)?shù)階微積分即非整數(shù)階微積分，數(shù)學(xué)家從不同的角度入手，給出了分?jǐn)?shù)階微積分不同的定義。每種定義有其表達(dá)式與性質(zhì)，其中被多數(shù)人所采納的定義形式是 Riemann-Liouville（R-L）分?jǐn)?shù)階微積分如下［11］：

式中：Γ為伽馬函數(shù)；n為任意整數(shù)；aDqt為分?jǐn)?shù)階微分算子；q為分?jǐn)?shù)階階次；a為積分算子下限；t為積分算子上限。

在工程中應(yīng)用較多及本文采用的是更為方便的Caputo定義，這種計(jì)算形式避免了極限求解，且初始條件以整數(shù)階微積分的形式給出，其定義如下：

（R－L）分?jǐn)?shù)階微積分定義與Caputo分?jǐn)?shù)階微積分定義都是在時域范圍內(nèi)對函數(shù)f（x）的定義。（RL）的Laplace變換涉及到分?jǐn)?shù)階積分的初始值與原（R－L）分?jǐn)?shù)微分的初始值，雖可求出含（R－L）分?jǐn)?shù)階次導(dǎo)數(shù)的解，但很難解賦予合理的物理意義。而Caputo分?jǐn)?shù)階微積分定義的最大優(yōu)點(diǎn)是其初始值與整數(shù)階微積分一樣，有明確的物理意義，具有超奇異性。

2．2 油氣懸架的分?jǐn)?shù)階模型

如前所述，影響懸架輸出特性的因素眾多，試驗(yàn)條件也無法做到規(guī)避所有干擾項(xiàng)，建立數(shù)學(xué)模型和模擬分析前，都需根據(jù)具體情況做出相應(yīng)假設(shè)。本文研究油氣懸架阻尼力等效分?jǐn)?shù)階黏彈性的表現(xiàn)，可借助筒式液力減振器的研究經(jīng)驗(yàn)［12］，對工作參數(shù)做了適當(dāng)?shù)暮喕秃雎?，做出如下假設(shè)：

（1）懸架中活塞、活塞桿與缸筒之間不產(chǎn)生泄漏；

（2）液壓油的重力遠(yuǎn)小于其工作情況下的運(yùn)動力，故不計(jì)工作油液的重力勢能；

（3）不計(jì)壓力變化而引起系統(tǒng)構(gòu)件的彈性變形；

（4）被研究封閉區(qū)域內(nèi)同一瞬時壓力處處相等；

（5）不考慮油氣懸架內(nèi)部溫度變化對性能的影響，即假設(shè)其工作工程中是恒溫過程。

油氣懸架的受力分析對象為活塞桿筒組件，活塞桿在內(nèi)部做往復(fù)運(yùn)動。油氣懸架存在壓縮運(yùn)動和拉伸運(yùn)動，這兩種運(yùn)動中流體的運(yùn)動方向是相反，因此在對應(yīng)的兩種運(yùn)動狀態(tài)下對活塞桿筒組件進(jìn)行力平衡分析。油氣懸架在受到外部激勵時，活塞桿及活塞組件和缸筒之間要產(chǎn)生相對運(yùn)動。從輸出力的角度出發(fā)整體來分析其振動特點(diǎn)。油氣懸架在受到外部激勵時，活塞桿及活塞組件和缸筒之間產(chǎn)生相對運(yùn)動時，其輸出作用力主要由氣體彈性力、油液阻尼力及活塞和缸筒之間的摩擦力組成。輸出合力公式如下：

式中：F為油氣懸架缸筒輸出力；Fk為氣體彈性力項(xiàng)；Fc為油液阻尼力項(xiàng)；Ff為活塞和缸筒之間密封件的摩擦力項(xiàng)

對于彈性力、阻尼力及摩擦力項(xiàng)的計(jì)算，先取基本整數(shù)階常用計(jì)算方法［13］進(jìn)行說明，彈性力項(xiàng)為：

式中：P為工作腔氣體壓力 ；A3為活塞桿面積；M為流體密度；γ為氣體多變指數(shù)；V0為靜平衡位置氣體體積；x為活塞和缸筒間相對位移。

油氣懸架密封件的摩擦力使得油氣懸架的動靜太特性呈現(xiàn)出遲滯性，物理特性類似阻尼可以借助液壓缸及活塞桿的密封件摩擦阻力的算法來計(jì)算：

式中：μ1為密封圈摩擦系數(shù)；ΔP為密封圈兩側(cè)壓力差；D為液壓缸缸筒內(nèi)徑d為活塞桿外徑；bD，bd為活塞、活塞桿密封圈寬度；kDkd為密封圈摩擦修正系數(shù)。

油氣懸架阻尼力項(xiàng)根據(jù)懸架阻尼孔的節(jié)流孔口形式，選取薄壁孔口類型的阻尼流量計(jì)算方法如下：

式中：Cz，Az為阻尼孔的流量系數(shù)、過流面積；Cd，Ad為工作腔氣體壓力；ρ為油液密度。

sign為符號函數(shù)，根據(jù)輸入信號的速度方向取值，壓縮行程為1，伸張行程－1。

先只考慮壓縮行程即sigx＝1。聯(lián)立方程（4）～（6）可以得到懸架輸出力整數(shù)階方程，便于計(jì)算進(jìn)行無量 綱 化 處 理，令 a1＝ （A1－A3）3r／［2n2（CzAz＋CdAd）2］，a2＝Mg·Vγ0，a3＝V0，a4＝A3，得到簡化系統(tǒng)方程：

根據(jù)圖1運(yùn)行原理及其幾何關(guān)系可知：

式中：z（t）為懸架上端輸出位移；y（t）為懸架下方輸入激勵位移。

上式描繪了三者位移及加速度關(guān)系，在油氣懸架的下方施加一個正弦位移激勵 y（t）＝A sinωt作為輸入，相當(dāng)于路面激勵，A為激勵振幅，ω為激勵頻率。

聯(lián)立式（3）～（9），在整數(shù)階方程基礎(chǔ)上對于耗散阻尼項(xiàng)引入分?jǐn)?shù)階參數(shù)，同時計(jì)入自身重力，建立系統(tǒng)等效分?jǐn)?shù)階 Bagley-Torvik阻尼振動方程［14］，進(jìn)一步得到油氣懸架輸出力的非線性分?jǐn)?shù)階方程：

式中：D為分?jǐn)?shù)階微分算子；α為分?jǐn)?shù)階算子階次。

方程（10）為分?jǐn)?shù)階非線性常微分方程，彈性力項(xiàng)中變量x位于分母上，阻尼力項(xiàng)亦為平方是非線性項(xiàng)，這使得懸架的振動輸出力是非線性的。大部分分?jǐn)?shù)階微分方程難以求得解析解，此類分?jǐn)?shù)階微分方程的解析解也只能由比較特殊的函數(shù)來表示，針對工程需求一般求取此類方程的數(shù)值解。

目前分?jǐn)?shù)階微積分常用的計(jì)算方法包括解析法、數(shù)值法和濾波器算法。本文采用改進(jìn)的Oustaloup濾波器［15］算法來仿真未知信號的分?jǐn)?shù)階微積分，其原理是用高階的傳遞函數(shù)近似替代分?jǐn)?shù)階算子。Oustaloup算法就是在頻率段內(nèi)，用泰勒展開并取一階近似得到傳遞函數(shù)。另一方面，考慮到濾波器分子和分母階次的一致性會導(dǎo)致仿真過程出現(xiàn)代數(shù)環(huán)，故而進(jìn)行改進(jìn)，后面再接一個低通濾波器。這樣可編寫濾波器零極點(diǎn)及增益函數(shù)，在Simulink中建立起如圖3的分?jǐn)?shù)階微分器模塊。進(jìn)而搭建懸架系統(tǒng)微分方程，求得數(shù)值解。

圖3 分?jǐn)?shù)階微分濾波器Fig．3 Fractional differential filter

3 懸架分?jǐn)?shù)階模型數(shù)值求解和試驗(yàn)

3．1 分?jǐn)?shù)階求解及試驗(yàn)對比

自行設(shè)計(jì)油氣懸架試驗(yàn)臺能獲取油氣懸架輸出力特性的動態(tài)數(shù)據(jù)，方便于對不同工況下的工作特性進(jìn)行分析和比較，驗(yàn)證油氣懸掛分?jǐn)?shù)階模型的正確性。文中試驗(yàn)使用的試驗(yàn)臺其結(jié)構(gòu)圖如4所示，主要由驅(qū)動裝置、懸架系統(tǒng)、傳感器、采集系統(tǒng)和簧上負(fù)載組成。設(shè)計(jì)思路主要是參照1／4車的振動模型，模擬車輛實(shí)際運(yùn)行過程中油氣懸架的工作狀態(tài)，測試其受到外部激勵和載荷時的輸出特性。

整個臺架驅(qū)動力由下端地基平臺上的變頻三相異步電動機(jī)發(fā)出，經(jīng)減速器后將力傳導(dǎo)至與支架接觸的偏心輪上。偏心輪與傳動梁保持接觸，傳動梁通過關(guān)節(jié)軸承與懸架活塞桿鉸接。偏心輪隨電機(jī)勻速旋轉(zhuǎn)同時驅(qū)動傳動梁做上下方向的擺動，擺動傳遞至懸架缸活塞桿與傳動梁鉸接處按照幾何關(guān)系即為豎直方向的正弦激勵。為防止試驗(yàn)時質(zhì)量塊左右晃動，在質(zhì)量塊兩側(cè)設(shè)計(jì)了限位支架。本設(shè)計(jì)激勵頻率可通過調(diào)節(jié)電機(jī)頻率來隨時改變；激勵幅值可通過調(diào)節(jié)偏心輪的偏心距來方便的調(diào)整；懸架簧上質(zhì)量可通過改變質(zhì)量塊數(shù)量來調(diào)整大小。

對于橫梁臂自身質(zhì)量及質(zhì)量塊質(zhì)量需通過杠桿比進(jìn)行等效，即為簧上質(zhì)量部分，通過氣室靜態(tài)氣壓推算簧上質(zhì)量予以校正。試驗(yàn)臺系統(tǒng)構(gòu)件引起的彈性變形以及系統(tǒng)自身各處阻尼較小可以忽略不計(jì)。試驗(yàn)臺下放為鋼結(jié)構(gòu)試驗(yàn)臺地基，地基自身振動引發(fā)的誤差極小也可忽略不計(jì)。每次試驗(yàn)前需對偏心距進(jìn)行校正，保證激勵振幅準(zhǔn)確性。傳動梁及橫臂梁的杠桿比問題，需首先參照幾何關(guān)系進(jìn)行等效轉(zhuǎn)換和試驗(yàn)校正。

圖4 油氣懸架試驗(yàn)臺Fig．4 Test stand of hydro-pneumatic suspension

為驗(yàn)證分?jǐn)?shù)階理論的準(zhǔn)確性與整數(shù)階對比，完整描述試驗(yàn)過程及結(jié)果，本文同時根據(jù)試驗(yàn)臺結(jié)構(gòu)設(shè)計(jì)構(gòu)建了臺架三維模型如圖5，進(jìn)而將模型導(dǎo)入多體動力學(xué)軟件ADAMS／View中，通過設(shè)置各零部件的材料屬性并添加零部件約束及力關(guān)系，與LMSAmesim液壓軟件聯(lián)立得到試驗(yàn)臺多剛體整數(shù)階仿真模型，目的是作為整數(shù)階結(jié)果與論文的分?jǐn)?shù)階模型進(jìn)行對比。

試驗(yàn)選用懸架基本參數(shù)如表1所示，為某噸位自卸車實(shí)車懸架。研究主要對比相同工況下懸架缸的輸出特性并針對不同工況展開試驗(yàn)，即簧上載荷的振動加速度、懸架缸輸出力幅值和頻率、運(yùn)動行程等。

圖5 試驗(yàn)臺架軟件仿真Fig．5 Software simulation of test stand

表1 試驗(yàn)用懸架基本參數(shù)Tab．1 The basic parameters of test suspension

本文針對偏心輪偏心距為5 mm，轉(zhuǎn)動頻率為1 Hz正弦激勵情況先進(jìn)行分析，這也是工程車輛常見激勵頻率。取多個不同分?jǐn)?shù)階次于振動穩(wěn)態(tài)后在時域范圍內(nèi)進(jìn)行對比分析發(fā)現(xiàn)，與整數(shù)階次結(jié)果的周期和相位基本一致，而幅值差異明顯，可見分?jǐn)?shù)階次的選取對振動結(jié)果的描述有很大影響。為了最為吻合試驗(yàn)數(shù)據(jù)，分析前述方程，將最優(yōu)階次的選取轉(zhuǎn)化為單變量的最優(yōu)數(shù)學(xué)問題［16］，單變量為分?jǐn)?shù)階次α，根據(jù)計(jì)算模型與實(shí)測數(shù)據(jù)間誤差最小原則進(jìn)行數(shù)值求解選取最優(yōu)階次。定義誤差評價函數(shù)E為計(jì)算結(jié)果與實(shí)測數(shù)據(jù)的差方均值，這樣更為直觀方便，取毫米為單位。如下：

式中：E為誤差評定；n為對比采樣點(diǎn)個數(shù)；xli為整數(shù)階結(jié)果第i個采樣點(diǎn)位移；xi為實(shí)測結(jié)果第i個采樣點(diǎn)位移。

同樣x2i為分?jǐn)?shù)階結(jié)果中第i個采樣點(diǎn)位移，代入后則為分?jǐn)?shù)階計(jì)算結(jié)果與實(shí)際測量誤差。約束條件的選取設(shè)在頻段［0．4，1．4］內(nèi)進(jìn)行，迭代步長為 0．001，觀察最優(yōu)階次變化后目標(biāo)函數(shù)的變化如圖7。迭代求解得分?jǐn)?shù)階建模在分?jǐn)?shù)階次取0．752～1的范圍內(nèi)比整數(shù)階結(jié)果誤差要小，在1處即為整數(shù)階結(jié)果同為1．68的差方均值。迭代結(jié)果最小誤差即最優(yōu)值出現(xiàn)在0．912處，其幅值結(jié)果最接近實(shí)測結(jié)果。下文以0．912次分?jǐn)?shù)階結(jié)果、整數(shù)階次仿真結(jié)果、實(shí)測數(shù)據(jù)進(jìn)一步對比分析。

圖6 差方均值隨迭代階次變化Fig．6 Error changes with the order

實(shí)測數(shù)據(jù)是懸架穩(wěn)定振動后測得，為了便于說明對比圖的繪制中，實(shí)測數(shù)據(jù)相位滯后或提前90°，取Matlab分?jǐn)?shù)階數(shù)值結(jié)果、ADAMS整數(shù)階結(jié)果及實(shí)測結(jié)果合并對比，如圖7為懸架缸缸筒相對位移對比，同取在4 s時域范圍內(nèi)對比，后續(xù)分析也如是。

分析圖7可知缸筒與活塞桿的相對位移中分?jǐn)?shù)階建模結(jié)果更加貼近于實(shí)際情況，整數(shù)階結(jié)果振幅比試驗(yàn)數(shù)據(jù)振幅要更大，最大振幅要超過實(shí)測值約1．5 mm。在分?jǐn)?shù)階原理上，對于阻尼項(xiàng)階數(shù)與運(yùn)動關(guān)系可知，隨著分?jǐn)?shù)階次的增大振動會越來越激烈，振動位移幅度也會越大，即前面已有的運(yùn)動狀態(tài)對下一時刻影響力減弱。反之越小于1，則增加這種相互關(guān)系，使得阻尼力變大，振幅減小。明顯分?jǐn)?shù)階建模在幅值方面更為精確，尤其是在位移下止點(diǎn)位置。受阻尼力影響，活塞壓縮和伸張行程是一個不重合過程。在伸張行程，在阻尼的作用下使得其伸張行程要更小，分?jǐn)?shù)階的引入強(qiáng)調(diào)了阻尼的非線性，分?jǐn)?shù)階模型結(jié)果在伸張后半程和壓縮前半程的優(yōu)勢更為明顯，最大值最小值的出現(xiàn)時間與實(shí)測結(jié)果更為相近，改變趨勢更為吻合。

圖7 缸筒相對位移對比Fig．7 data of displacement curve comparison diagram

圖8 缸筒加速度對比Fig．8 Accelerated speed of cylinder comparison

圖9 懸架缸輸出力曲線對比Fig．9 Output force of suspension cylinders comparison

圖8 為懸架缸缸筒加速度曲線，與圖6相似，頻率和幅值都很接近，明顯仿真結(jié)果的曲線比試驗(yàn)結(jié)果更為平滑，這是因?yàn)樵谡麛?shù)階仿真模型及軟件中把試驗(yàn)臺當(dāng)作剛體處理，沒有考慮傳動梁等構(gòu)件自身的阻尼及多處振動，但其對結(jié)果影響不大。

圖9為懸架缸輸出力曲線，動平衡后油氣懸架的輸出力呈周期性變化，結(jié)果曲線幅值變化頻率與實(shí)測數(shù)據(jù)相近。油氣懸架本身的高度耦合多復(fù)雜工況使得黏彈性特性出現(xiàn)且最終趨向穩(wěn)態(tài)。在壓縮行程分?jǐn)?shù)階模型優(yōu)勢并不明顯，而在伸張行程分?jǐn)?shù)階建模遠(yuǎn)比整數(shù)階結(jié)果更貼近實(shí)測數(shù)據(jù)，在伸張行程主要是阻尼力在起作用，分?jǐn)?shù)階模型的引入是針對阻尼力進(jìn)行了調(diào)整，其優(yōu)勢在伸張行程體現(xiàn)明顯，阻尼力的調(diào)整同樣使得輸出力較之于整數(shù)階相位上有延遲，其輸出力最大值出現(xiàn)要晚，這與實(shí)測數(shù)據(jù)更為相近。

3．2 多工況下分?jǐn)?shù)階次與誤差變化

為了進(jìn)一步驗(yàn)證分?jǐn)?shù)階建模在油氣懸架多工況中具有廣泛的有效性，可改變激勵振幅和頻率，觀察最優(yōu)分?jǐn)?shù)階次的變化，同時對比其與整數(shù)階建模誤差變化。本文通過控制偏心輪距和激勵電機(jī)的轉(zhuǎn)速改變激勵頻率和振幅，進(jìn)行多次試驗(yàn)，并參照分?jǐn)?shù)階結(jié)果差方均值最小為目標(biāo)進(jìn)行最優(yōu)階次的確定。在同一振幅5 mm改變激勵頻率，在同一頻率1 Hz下改變振幅分別進(jìn)行最優(yōu)化求解和試驗(yàn)，通過定量計(jì)算擬合曲線得到圖10最優(yōu)階次隨激勵頻率變化以及圖11隨激勵振幅變化。

分析圖10和11發(fā)現(xiàn)隨著頻率的改變最優(yōu)階次下降并不大，高頻激勵下最終穩(wěn)定趨于0．9，而隨著激勵振幅的變化下降較大趨于穩(wěn)定0．86，這是因?yàn)樵诟叻焖僬駝酉伦枘崃ζ鸬搅烁鼮槊黠@的作用，其等效黏彈性明顯，二者最后都趨于穩(wěn)定，這說明分?jǐn)?shù)階建模的取值在高幅階段可以穩(wěn)定選取，建模方便直接。

圖10 最優(yōu)階次隨激勵頻率變化Fig．10 Optimal order changes with the excitation frequency

圖11 最優(yōu)階次隨激勵振幅變化Fig．11 Optimal order changes with the excitation amplitude

圖12 不同振幅下差方均值對比Fig．12 Error comparison at different amplitude

為了進(jìn)一步分析和觀察分?jǐn)?shù)階最優(yōu)階次確定后與整數(shù)階建模的差別，可對二者按評價指標(biāo)進(jìn)行觀察。取多組不同振幅下數(shù)據(jù)進(jìn)行評價然后擬合繪圖得到圖12即不同振幅激勵下差方均值變化，差方的平均值越小說明精確度越高，可見隨著振幅增大二者誤差都有進(jìn)一步的增大，但是整數(shù)階誤差增大更為明顯，進(jìn)一步驗(yàn)證了分?jǐn)?shù)階建模的有效性。另一方面，取多組不同頻率下試驗(yàn)數(shù)據(jù)擬合可得不同頻率下差方均值對比，因最優(yōu)階次下降不大，二者誤差并不明顯，但在總體范圍內(nèi)變頻和變幅值的分?jǐn)?shù)階建模誤差都要小于整數(shù)階。結(jié)合圖11可知等效黏彈性的分?jǐn)?shù)階建模方式能進(jìn)行有效修正，得到更準(zhǔn)確的結(jié)果。

圖13 阻尼力對比Fig．13 Damping force comparison

分析分?jǐn)?shù)階建模的特點(diǎn)可以對阻尼力項(xiàng)進(jìn)行說明，圖13為初始試驗(yàn)條件即1 Hz和5 mm激勵下二者數(shù)值計(jì)算的阻尼力項(xiàng)對比，在伸張行程分?jǐn)?shù)階可以達(dá)到穩(wěn)定的550 N，而整數(shù)階結(jié)果只有400 N，分?jǐn)?shù)階阻尼對系統(tǒng)消耗更大，而二者彈性力項(xiàng)變化不大，原理上這造成了懸架的相對位移非標(biāo)準(zhǔn)正弦波，在阻尼力的非線性作用下使得向下運(yùn)動的伸張行程偏大，向上運(yùn)動的壓縮行程偏小。

4 結(jié) 論

（1）根據(jù)油氣懸架工作原理，突出油氣懸架的多相介質(zhì)特點(diǎn)，采用等效黏彈性阻尼的分析方法，引入分?jǐn)?shù)階導(dǎo)數(shù)建立了油氣懸架的 Bagley-Torvik方程。利用Oustaloup算法對分?jǐn)?shù)階非線性微分方程進(jìn)行數(shù)值求解得到油氣懸架分?jǐn)?shù)階輸出力，位移，加速度等多項(xiàng)數(shù)值結(jié)果。

（2）將分?jǐn)?shù)階、整數(shù)階仿真與試驗(yàn)結(jié)果進(jìn)行對比，驗(yàn)證了分?jǐn)?shù)階模型的正確性。在1 Hz，5 mm激勵下測試表明取分?jǐn)?shù)階次0．912時與實(shí)驗(yàn)結(jié)果最為相近，該簡潔模型能較精確地描述油氣懸架的特性。

（3）改變激勵振幅和頻率，在不同工況下進(jìn)行對比分析最優(yōu)階次的變化，在高頻激勵下最優(yōu)階次趨于0．9，在高幅激勵下趨于0．86，在整個測試激勵范圍內(nèi)分?jǐn)?shù)階模型誤差都要小于整數(shù)階，這說明了分?jǐn)?shù)階建模的有效性。

［1］Pazooki A，Cao Dong-pu，Rakheja S．Experimental and analytical evaluations of a torsio-elastic suspension for off-road vehicles［J］．SAE International Journal of Materials and Manufacturing，2010，3（1）：326－338．

［2］Gao B，Darling J，Tilley D，et al．Control of a hydropneumatic active suspension based on a non－linear quarter-car model［J］．Journal of Systems＆Control Engineering，2006，220（1）：15－31．

［3］孫濤，喻凡，張振東．重載車輛油氣懸架不確定分析及H∞控制器設(shè)計(jì)［J］．振動與沖擊，2007，26（9）：51－54．SUN Tao，YU Fan，ZHANGZhen-dong．Uncertainty analysis for a hydro-pneumatic suspension system of heavy duty vehicles and H∞controller design［J］．Journal of Vibration and Shock，2007，26（9）：51－54．

［4］楊波，陳思忠，王勛，等．雙氣室油氣懸架特性研究［J］．機(jī)械工程學(xué)報，2009，45（5）：276－280．YANG Bo，CHEN Sizhong，WANG Xun，et al．Research of Twin accumulator Hydro－pneumatic Suspension［J］．Journal of Mechanical Engineerin，2009，45（5）：276－280．

［5］王增全，申焱華，楊玨．連通式油氣懸架數(shù)學(xué)模型及特性分析［J］．農(nóng)業(yè)工程學(xué)報，2012，28（5）：60－65．WANG Zeng-quan，SHEN Yan-hua，YANG Jue．Mathematical model and characteristics analysis of interconnected hydropneumatic suspension［J］．Transactions of the Chinese Society of Agricultural Engineering，2012，28（5）：60－65．

［6］徐道臨，張林，周加喜．重型礦用自卸車油氣懸架參數(shù)優(yōu)化［J］．振動與沖擊，2012，31（24）：98－101．XU Dao-lin， ZHANG Lin， ZHOU Jia-xi． Parametric optimization of hydro pneumatic suspension of a heavy mining dumper［J］．Journal of Vibration and Shock，2012，31（24）：98－101．

［7］黃夏旭，申焱華，楊玨，等．基于集中參數(shù)熱模型法的自卸車油氣懸架系統(tǒng)熱分析［J］．農(nóng)業(yè)工程學(xué)報，2013，29（10）：64－70．HUANG Xia-xu， SHEN Yan-hua， YANG Jue， et al．Thermal analysis of hydro-pneumatic suspension system for dumper based on a lumped-parameter thermal model［J］．Transactions of the Chinese Society of Agricultural Engineering，2013，29（10）：64－70．

［8］張輝，何景峰，韓俊偉．油氣懸架加載試驗(yàn)臺正弦波復(fù)現(xiàn)控制方法［J］．農(nóng)業(yè)機(jī)械學(xué)報，2009，40（7）：35－39．ZHANG Hui，HE Jing-feng， HAN Jun-wei． Sine wave replication control methodon hydro pneumatic suspension force loading system［J］．Transactions of the Chinese Society for Agricultural Machinery，2009，40（7）：35－39．

［9］黃夏旭，楊玨，申焱華，等．基于氣體溶解與油液可壓縮性的油氣懸架性能研究［J］．農(nóng)業(yè)機(jī)械學(xué)報，2013，44（6）：14－18．HANG Xia-xu， YANG Jue， SHEN Yan-hua， et al．Characteristics analysis of hydro pneumatic suspension based on gas dissolution and oil compressibility［J］．Transactions of the Chinese Society for Agricultural Machinery， 2013，44（6）：14－18．

［10］ Song D Y，J TQ．Study on the constitutive equation with fractonal derivative for viscoelastic fluids－modified Jeffreys models and its application［J］．Rhelo Acta 1998，21（37）：512－517．

［11］Connolly J A．The numerical solution of fractional and distributed order differential equations［D］．University of LIverpool，2004．

［12］ Kurfess，Thomas R，Messman，et al．A methodology to integrate a nonlinear shock absorber dynamics into a vehicle model for system identification［J］．SAE International Journal of Materials and Manufacturing，2011，4（1）：527－534．

［13］ Wilkinson，Peter M，Van Dierendonck，et al．Theoretical model for the influence of gas properties and pressure on single-bubble formation at an orifice［J］．Source：Chemical Engineering Science，1994，9（49）：1429－1438．

［14］沈淑君，劉發(fā)旺．解分?jǐn)?shù)階Bagley-Torvik方程的一種計(jì)算有效的數(shù)值方法［J］．廈門大學(xué)學(xué)報（自然科學(xué)版），2004，7（3）：306－311．SHEN Shu-jun， LIU Fa-wang． A computational efficient numerical method of solving fractional Bagley－Torvik equation［J］．Journal of Xiamen University（Natural Science）2004，7（3）：306－311．

［15］薛定宇．高等應(yīng)用數(shù)學(xué)問題的MATLAB求解［M］．北京：清華大學(xué)出版社，2008：398－399．

［16］劉鯖潔，陳桂明，楊旗．基于Matlab工具的遺傳算法求解有約束最優(yōu)化問題［J］．兵工自動化，2008，27（11）：43－44．LIU Qing-jie， CHEN Gui-ming， YANG Qi． Solving constrained optimization through genetic algorithm based on matlab toolbox ［J］． Advanced Manufacture and Management，2008，27（11）：43－44．