農(nóng)用車輛橡膠隔振元件動力學(xué)模型研究*

高琦，邱緒云，宋裕民

(山東交通學(xué)院汽車工程學(xué)院，濟(jì)南市，250357)

通訊作者：邱緒云，男，1977年生，山東臨沂人，博士，教授；研究方向為車輛底盤設(shè)計與控制。E-mail: qiuxuyun@163.com

0 引言

我國農(nóng)業(yè)機(jī)械化迅猛發(fā)展的過程中，農(nóng)用車輛的廣泛應(yīng)用大大提高了農(nóng)業(yè)生產(chǎn)和運(yùn)輸?shù)男蔥1]。農(nóng)用車輛的使用工況相對惡劣，在作業(yè)和行駛過程中使車輛各結(jié)構(gòu)部件產(chǎn)生劇烈振動，同時發(fā)動機(jī)產(chǎn)生的振動也會傳遞給車輛各系統(tǒng)及其部件，不但會降低農(nóng)用車輛作業(yè)穩(wěn)定性，引起結(jié)構(gòu)噪聲，也會影響各部件的疲勞壽命[2-3]。

橡膠隔振元件作為農(nóng)用車輛動力裝置、減振裝置等部件間的連接元件與隔振元件，其動力學(xué)特性直接影響農(nóng)用車輛的隔振與降噪性能，它具有很強(qiáng)的非線性粘彈力學(xué)特性，受外界因素(加載幅值、加載頻率、工作周期以及溫度等)的影響比較大[4]。橡膠隔振元件模型的精度是農(nóng)用車輛各子系統(tǒng)及整車動力學(xué)仿真精度的關(guān)鍵影響因素之一，目前尚沒有公認(rèn)可接受的高精度橡膠隔振元件模型用于整車及子系統(tǒng)模型動態(tài)仿真。因此，建立能夠準(zhǔn)確描述橡膠隔振元件動態(tài)特性的動力學(xué)模型，對于提高農(nóng)用車輛整車模型仿真精度以及準(zhǔn)確預(yù)測復(fù)雜工況下各子系統(tǒng)特性對車輛隔振與降噪性能的影響具有重要意義。

目前為止，國內(nèi)外學(xué)者提出了不同的動力學(xué)模型用于描述橡膠隔振元件的動態(tài)力學(xué)特性，為克服Kelvin-Voigt模型、Maxwell模型、三參數(shù)Maxwell模型、廣義Maxwell模型以及Zener模型的不足，分?jǐn)?shù)導(dǎo)數(shù)模型因模型參數(shù)相對較少、調(diào)整能力強(qiáng)、能反映較寬頻域內(nèi)加載歷程對橡膠材料動態(tài)特性的影響而日益受到關(guān)注[5-10]。這類模型多為階次不大于1的低階分?jǐn)?shù)導(dǎo)數(shù)模型，在預(yù)測橡膠材料的能量損耗頻率相關(guān)性方面存在一定誤差，建立的模型往往需要進(jìn)行一定的修正和改進(jìn)才能更加精確地試驗結(jié)果吻合。最高階次大于1的高階分?jǐn)?shù)導(dǎo)數(shù)模型更接近于描述橡膠材料物理本質(zhì)[11-13]，但相關(guān)研究相對較少，因此，本文基于高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性理論建立彈性單元、摩擦單元和高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性單元疊加的橡膠隔振元件動力學(xué)模型，提出基于廣義多維學(xué)習(xí)策略和自適應(yīng)調(diào)整策略的粒子群算法的粘彈性單元參數(shù)辨識方法，并推導(dǎo)高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性單元的數(shù)值求解方法，為實現(xiàn)橡膠隔振元件在農(nóng)用車輛各系統(tǒng)及整車動力學(xué)仿真分析中的精確表征奠定基礎(chǔ)。

1 橡膠隔振元件動態(tài)力學(xué)試驗

選取某農(nóng)用車輛發(fā)動機(jī)前懸置軟墊為試驗對象，基于MTS831襯套試驗機(jī)進(jìn)行了軸向動態(tài)加載試驗，如圖1所示。

圖1 某農(nóng)用車輛發(fā)動機(jī)前懸置軟墊軸向動態(tài)加載試驗Fig. 1 Axial dynamic loading test for the engine front mounting cushion of an agricultural vehicle

通過分析該橡膠軟墊的工作特性，先對其施加1 000 N 的預(yù)載，然后進(jìn)行頻率為1～50 Hz，振幅分別為0.1、0.2、0.4、1.0和2.0 mm的動態(tài)力學(xué)加載試驗，得到在不同加載振幅下其動剛度與阻尼系數(shù)隨加載頻率的變化曲線，如圖2和圖3所示。

圖2 不同加載幅值下橡膠軟墊動剛度隨頻率的變化Fig. 2 Variation of dynamic stiffness of rubber cushion with frequency under different loading amplitudes

圖3 不同加載幅值下橡膠軟墊阻尼系數(shù)隨頻率的變化Fig. 3 Variation of damping coefficient of rubber cushion with frequency under different loading amplitudes

從圖2和圖3可以看出，在1～50 Hz頻率范圍內(nèi)，在同一振幅下，橡膠軟墊動態(tài)剛度與阻尼系數(shù)均隨頻率的增加而增大；在頻率相同時，隨著振幅的增加，橡膠軟墊的動剛度逐漸下降，阻尼系數(shù)具有先增大再減小的趨勢。

2 橡膠隔振元件動力學(xué)模型

所建立的橡膠隔振元件動力學(xué)模型由彈性單元、摩擦單元以及粘彈性單元疊加組成，如圖4所示。Fe、Ff、Fv和F依次為橡膠隔振元件模型的彈性力、摩擦力、粘彈力和總響應(yīng)力，N。

圖4 橡膠隔振元件動力學(xué)模型示意圖Fig. 4 Dynamic model of rubber vibration isolation component

2.1 彈性單元

彈性單元用來描述橡膠隔振元件的靜態(tài)主剛度。橡膠隔振元件的彈性特性在線變形范圍內(nèi)符合胡克定律，當(dāng)變形量超過一定范圍后，其彈性特性呈現(xiàn)出明顯的非線性。彈性單元中，彈性力與加載位移的關(guān)系可使用彈簧模型來描述。文獻(xiàn)[14]指出，與正切彈簧模型和結(jié)構(gòu)彈性單元模型相比，冪次多項式模型可以靈活而有效地描述彈性力與加載位移之間的線性或非線性關(guān)系。因此，選用冪次多項式模型描述橡膠隔振元件的彈性特性。

Fe=a0+a1x+a2x2+a3x3+…+anxn

(1)

式中：Fe——彈性力，N；

x——加載位移，mm；

a0,a1,a2,a3,…,an——冪次項的系數(shù)。

在振幅為x0的簡諧激勵下，彈性力的幅值

(2)

2.2 摩擦單元

摩擦單元中，摩擦力與加載位移的關(guān)系[13]

(3)

式中：Ff——摩擦力，N；

Ffmax——最大摩擦力，N；

x2——Ff從0開始，逐漸增大至Ffmax/2時的位移，mm；

xs——摩擦力與位移變化曲線上的起始點(diǎn)或反向點(diǎn)的位移，mm；

Ffs——摩擦力與位移變化曲線上位移為xs的摩擦力，N；

(xs，F(xiàn)fs)——摩擦力與位移變化曲線上的狀態(tài)參考點(diǎn)，此處取(0，0)。

在振幅為x0的簡諧激勵下，摩擦力的幅值

(4)

一個加載循環(huán)損耗的能量為

(5)

其中u0=Ff0/Ffmax。

2.3 粘彈性單元

選用分?jǐn)?shù)階Kelvin-Voigt模型與分?jǐn)?shù)階Maxwell模型并聯(lián)的高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性模型描述橡膠材料的粘彈性特性[13]，如圖5所示。圖5中k1和c1分別為分?jǐn)?shù)階Kelvin-Voigt模型中彈簧的彈性模量和黏壺的粘性系數(shù)，k2和c2分別為分?jǐn)?shù)階Maxwell模型中彈簧的彈性模量和黏壺的粘性系數(shù)。k1和k2用于反映材料的彈性模量，c1和c2用于反映材料的粘性。

圖5 高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性模型示意圖Fig. 5 High-order fractional derivative viscoelastic model

根據(jù)圖5所示模型，可得粘彈力與位移的關(guān)系

λ1λ2x(t)

(6)

其中λ1=k1/c1，λ2=k2/c2。

式中：Fv(t)——粘彈力，N；

x(t)——加載位移，mm；

α、β、γ——取值范圍(0，1)的分?jǐn)?shù)導(dǎo)數(shù)階數(shù)，其中α+γ的最高階數(shù)為2。

對式(6)進(jìn)行傅里葉變換，得該分?jǐn)?shù)導(dǎo)數(shù)粘彈力模型的復(fù)剛度

(7)

(8)

(9)

一個加載循環(huán)損耗的能量為

(10)

2.4 橡膠隔振元件動力學(xué)模型

將彈性單元、摩擦單元與分?jǐn)?shù)導(dǎo)數(shù)粘彈性單元疊加，得到橡膠隔振元件動力學(xué)模型的總響應(yīng)力與加載位移的關(guān)系

F(x)=Fe(x)+Ff(x)+Fv(x)

(11)

式中：F(x)——橡膠隔振元件動力學(xué)模型的總響應(yīng)力，N；

Fe(x)——彈性力，N；

Ff(x)——摩擦力，N；

Fv(x)——粘彈力，N。

在x=x0sin(ωt)的簡諧激勵下，總響應(yīng)力振幅

(12)

橡膠隔振元件模型的動剛度

(13)

橡膠隔振元件模型的阻尼系數(shù)

(14)

3 橡膠隔振元件動力學(xué)模型參數(shù)辨識

由于建立的橡膠隔振元件動力學(xué)模型，各單元之間滿足疊加原理，因此，模型參數(shù)的獲取可依次通過辨識摩擦單元、彈性單元以及粘彈性單元的參數(shù)實現(xiàn)。

3.1 摩擦單元參數(shù)辨識

由于加載速度較低時，粘彈性單元對橡膠襯套的響應(yīng)特性影響較小，為保證摩擦力能夠得到充分體現(xiàn)，采用低頻率、大振幅的動態(tài)加載工況對摩擦單元參數(shù)進(jìn)行辨識，摩擦單元參數(shù)辨識示意圖如圖6所示。

圖6 低頻加載工況下摩擦單元參數(shù)辨識示意圖Fig. 6 Parameter identification of friction element under low frequency loading

圖6中，加載位移達(dá)到最大時，摩擦力穩(wěn)定達(dá)到最大，此時曲線上端面接近極限位置處的斜率可近似表示該元件在此段位移的彈性剛度Ke。延長接近極限位置的上、下兩條切線，摩擦力模型中最大摩擦力Ffmax的值即為兩條切線間的垂直距離的1/2。取曲線下端面最大斜率的值Kmax，由式(15)可求得，摩擦單元中的x2[15]。

(15)

3.2 彈性單元參數(shù)辨識

彈性分力屬于橡膠隔振元件的靜態(tài)特性，也可通過低頻率、大振幅的動態(tài)加載工況對彈性單元的參數(shù)進(jìn)行辨識。為去除摩擦分力的影響，辨識出摩擦單元參數(shù)后，除去遲滯環(huán)曲線中摩擦分力分布，可得到彈性單元的剛度曲線，使用最小二乘法對彈性單元中多項式的系數(shù)進(jìn)行擬合，可得多項次冪和冪次項的系數(shù)。

3.3 粘彈性單元參數(shù)的辨識

橡膠隔振元件動態(tài)特性的頻率相關(guān)性主要通過粘彈性單元體現(xiàn)，為盡量減小摩擦力的影響，需通過小振幅加載工況獲取。為充分體現(xiàn)其在不同頻率工況下的粘彈性特性，選取振幅為0.1 mm，頻率別1、2、4、8、12、16、20、30、40和50 Hz共10個正弦加載試驗工況進(jìn)行粘彈性單元k1，c1，k2，c2，α，β和γ共7個參數(shù)的辨識。

粘彈性單元參數(shù)辨識實質(zhì)上是一個最優(yōu)參數(shù)估計問題，為使識別結(jié)果盡可能接近試驗測量結(jié)果，通過建立橡膠元件模型計算的動態(tài)響應(yīng)與試驗得到元件的動態(tài)響應(yīng)誤差最小的優(yōu)化函數(shù)來實現(xiàn)。此處基于橡膠隔振元件動剛度和阻尼系數(shù)誤差建立優(yōu)化目標(biāo)函數(shù)

(16)

式中：M——辨識工況的總個數(shù)；

由于粘彈性單元中包含高階分?jǐn)?shù)導(dǎo)數(shù)項，具有復(fù)雜的非線性，且需要辨識的參數(shù)較多，傳統(tǒng)的數(shù)學(xué)尋優(yōu)算法對非線性系統(tǒng)參數(shù)進(jìn)行辨識時，要求目標(biāo)函數(shù)連續(xù)可導(dǎo)，且搜尋最優(yōu)解時易陷入局部最優(yōu)。粒子群優(yōu)化(Particle Swarm Optimization，PSO)算法作為一種非解析仿生智能尋優(yōu)算法，對非線性系統(tǒng)參數(shù)辨識具有較好的應(yīng)用效果。為了提升傳統(tǒng)的PSO算法收斂速度較慢，降低陷入局部極值點(diǎn)的可能性，本文選取一種基于廣義多維學(xué)習(xí)策略和自適應(yīng)調(diào)整策略的PSO算法[16]對粘彈性單元的參數(shù)進(jìn)行辨識。

以式(16)所示目標(biāo)函數(shù)作為適應(yīng)度函數(shù)，算法中更新粒子的方法如式(17)所示。

i=1,2,3,…,n;j=1,2,3,…,d

(17)

式中：t——當(dāng)前迭代次數(shù)；

zi,j(t)——第t代中粒子適應(yīng)度值遞減順序排序后第i(1≤i≤n)個粒子Pi的位置向量的第j(1≤j≤d)維；

vi,j(t)——粒子Pi的速度向量的第j維；

bi(t)——一個隨機(jī)產(chǎn)生的概率；

(18)

式中：d——粒子的搜索維數(shù)；

n——群中的粒子數(shù)；

粒子Pi根據(jù)式(19)更新自己的每一維速度

(19)

式中：ωEt——動態(tài)調(diào)整的慣性權(quán)重；

zg,j(t)——粒子Pi的位置向量第j維的學(xué)習(xí)對象Pg的位置向量第j維；

κ——社會影響系數(shù)，κ取值比較大時，可能導(dǎo)致過早收斂到平均行為的現(xiàn)象；

r1、r2——廣義多維學(xué)習(xí)部分和社會學(xué)習(xí)部分的權(quán)重，取值為0～1之間均勻分布的隨機(jī)數(shù)。

(20)

式中：ni——落入每個小區(qū)間的粒子數(shù)。

為保證粘彈性單元參數(shù)識別的準(zhǔn)確性，提高模型的預(yù)測精度，參數(shù)識別過程中，需使模型計算的每個頻率工況下動剛度和阻尼系數(shù)，與試驗值均不能有較大誤差。通常認(rèn)為模型的計算精度在90%以上，可滿足工程應(yīng)用要求，為此，建立約束條件

(21)

參數(shù)識別過程中算法的收斂性準(zhǔn)則為

(22)

式中：σfobj——所有粒子適應(yīng)值的標(biāo)準(zhǔn)差；

μfobj——所有粒子適應(yīng)值的平均值。

該判據(jù)可以避免在fobj的全局最優(yōu)值附近出現(xiàn)小的波動。否則，當(dāng)?shù)螖?shù)達(dá)到最大迭代次數(shù)Gmax時，優(yōu)化過程終止。

經(jīng)多次調(diào)試，算法各參數(shù)設(shè)置如下：粒子維數(shù)d為7，位置向量x=(k1,k2,c1,c2,α,β,γ)，各維搜索空間下限zmin=(0,0,0,0,0,0,0)，上限zmax=(103,103,103,103,1,1,1)；群中的粒子數(shù)n取30，廣義多維學(xué)習(xí)部分權(quán)重r1和社會學(xué)習(xí)部分權(quán)重r2為0～1之間均勻分布的隨機(jī)數(shù)，最大迭代次數(shù)Gmax取100。粘彈性單元參數(shù)辨識的流程如圖7所示。

圖7 粘彈性單元參數(shù)辨識流程圖Fig. 7 Flow chart of identification for viscoelastic element parameters

3.4 橡膠隔振元件模型參數(shù)辨識結(jié)果與驗證

橡膠隔振元件模型各單元參數(shù)辨識結(jié)果如表1所示。各單元參數(shù)辨識完成后，基于橡膠隔振元件動力學(xué)模型進(jìn)行了正弦激勵動態(tài)仿真計算。將計算結(jié)果與試驗結(jié)果進(jìn)行對比，如圖8和圖9所示。

從圖8和圖9可以看出，模型仿真計算的橡膠隔振件的動剛度與阻尼系數(shù)的頻率響應(yīng)特性與幅值響應(yīng)特性均與試驗結(jié)果基本一致。由于摩擦單元損耗的能量不受頻率影響，使在激勵頻率較小時，模型仿真計算出的橡膠隔振能量損耗不為零。由此可見，所建模型能夠準(zhǔn)確地描述橡膠隔振元件動態(tài)特性的分布特征。

在振幅0.1～2 mm、頻率1～50 Hz范圍內(nèi)，對模型仿真計算結(jié)果與試驗結(jié)果進(jìn)行統(tǒng)計，得到該模型在不同振幅工況下，動剛度的平均誤差和最大誤差分別為3.08%和4.52%，阻尼系數(shù)的平均誤差和最大誤差分別為4.58%和5.79%，說明建立的模型在不同振幅和頻率下均有較高的精度。

表1 橡膠隔振元件模型各單元參數(shù)辨識結(jié)果Tab. 1 Identification results of element parameters of

圖8 橡膠隔振元件動剛度模型計算值與試驗值對比Fig. 8 Comparison of dynamic stiffness between model calculation value and test value of rubber vibration isolation element

圖9 橡膠隔振元件阻尼系數(shù)模型計算值與試驗值對比Fig. 9 Comparison of damping coefficients between model calculation value and test value of rubber vibration isolation element

4 橡膠隔振元件動力學(xué)模型數(shù)值求解

為了將建立的橡膠隔振元件動力學(xué)模型應(yīng)用于農(nóng)用車輛動力學(xué)建模與分析中，對橡膠隔振元件模型的數(shù)值求解方法進(jìn)行推導(dǎo)。

4.1 高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性單元數(shù)值求解

根據(jù)整數(shù)階微積分理論，可導(dǎo)函數(shù)f(t)對時間t的整數(shù)階導(dǎo)數(shù)可由各階向左做差商的極限定義，當(dāng)m取自然數(shù)時，可推導(dǎo)f(t)的m階導(dǎo)數(shù)為

(23)

式中： Δt——時間變量t在取值范圍內(nèi)微小增量。

將式(23)整數(shù)階導(dǎo)數(shù)的定義推廣到任意階，即將式(23)中的整數(shù)m用任意實數(shù)ε取代，并引入Gamma函數(shù)，可得到Grunwald-Letnikov形式的分?jǐn)?shù)階導(dǎo)數(shù)

f(t-jΔt)，ε>0

(24)

式中： Γ(·)——Gamma函數(shù)，Γ(n)=(n-1)!。

設(shè)函數(shù)f(t)定義在區(qū)間[0，T](T為大于0的任意實數(shù))上，用分?jǐn)?shù)代替Δt，即Δt=t/N(N為自然數(shù))，則式(24)可以寫成

(25)

由于

(26)

因此，可得到

(27)

零初始條件下，當(dāng)N取一極大值時，可進(jìn)一步得到分?jǐn)?shù)階導(dǎo)數(shù)的數(shù)值計算公式

(28)

(29)

由式(6)可知，橡膠隔振元件模型中粘彈性單元的微分方程含多個分?jǐn)?shù)導(dǎo)數(shù)項。一個微分方程中含有多個導(dǎo)數(shù)項的方程進(jìn)行數(shù)值計算，根據(jù)分?jǐn)?shù)導(dǎo)數(shù)運(yùn)算的可加性，可先把方程降階為多個方程組成的方程組，使每個方程中含有一個導(dǎo)數(shù)項，再進(jìn)行分?jǐn)?shù)階計算[17]。

考慮如下形式的分?jǐn)?shù)階微分方程

0

(30)

式中：x(t)——位移，mm；

Fs(t)——響應(yīng)力，N；

sj——分?jǐn)?shù)階導(dǎo)數(shù)項的系數(shù)；

sl——一次項的系數(shù)；

σl、σl-1——分?jǐn)?shù)階導(dǎo)數(shù)的階數(shù)，σl和σl-1為大于0的實數(shù)，且σl>σl-1。

式中：θh——取值在(0，1]區(qū)間內(nèi)的實數(shù)。

根據(jù)分?jǐn)?shù)導(dǎo)數(shù)運(yùn)算的可加性，對式(30)進(jìn)行降階處理，設(shè)

(31)

用y1替換x(t)，則式(31)可變形為

(32)

(33)

以橡膠隔振元件的加載位移作為輸入時，可先將位移的各導(dǎo)數(shù)項按上述方法進(jìn)行降階處理為式(32)所示的形式，并根據(jù)迭代公式(29)對方程組各項進(jìn)行數(shù)值計算，得到各位移項當(dāng)前時刻的值；然后再將粘彈力的各項按上述方法進(jìn)行降階處理，再結(jié)合迭代公式(29)，求解出當(dāng)前時刻粘彈力值。

4.2 數(shù)值求解結(jié)果驗證

在每一時刻的加載位移下，根據(jù)4.1節(jié)推導(dǎo)的數(shù)值求解方法求取粘彈力，再根據(jù)式(1)與式(3)分別求取彈性力與摩擦力，最后將上述各單元力代入式(11)，求得橡膠隔振元件模型的總響應(yīng)力。

分別取頻率1 Hz、振幅2 mm和頻率50 Hz、振幅1 mm兩種正弦信號作為位移輸入，進(jìn)行數(shù)值求解得到橡膠隔振元件的總響應(yīng)力，并與試驗結(jié)果進(jìn)行對比，如圖10和圖11所示。

圖10 1 Hz、2 mm工況下仿真與試驗力—位移曲線Fig. 10 Simulation and test force-displacement curves under 1 Hz and 2 mm working condition

圖11 50 Hz、1 mm工況下仿真與試驗力—位移曲線Fig. 11 Simulation and test force-displacement curves under 50 Hz and 1 mm working condition

從圖10和圖11可以看出，在不同頻率工況下，數(shù)值求解得到橡膠隔振元件的力—位移曲線與試驗測得的力—位移曲線具有較好的一致性。

5 結(jié)論

1) 以某農(nóng)用車輛發(fā)動機(jī)前懸置軟墊為載體，通過動態(tài)加載試驗對其軸向頻率響應(yīng)特性和幅值響應(yīng)特性進(jìn)行了研究，建立了彈性單元、摩擦單元和高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性單元疊加的橡膠隔振元件動力學(xué)模型。

2) 運(yùn)用基于廣義多維學(xué)習(xí)策略和自適應(yīng)調(diào)整策略的粒子群算法進(jìn)行了粘彈性單元參數(shù)辨識，結(jié)合試驗結(jié)果，辨識得到了模型參數(shù)。經(jīng)試驗驗證，在振幅0.1～2 mm、頻率1～50 Hz范圍內(nèi)，所建模型仿真計算與試驗得到動剛度的平均誤差和最大誤差分別為3.08%和4.52%，阻尼系數(shù)的平均誤差和最大誤差分別為4.58%和5.79%，說明建立的模型能夠精確地預(yù)測橡膠隔振元件的動態(tài)響應(yīng)特性，以期為復(fù)雜工況下準(zhǔn)確預(yù)測農(nóng)用車輛系統(tǒng)動態(tài)響應(yīng)特性奠定基礎(chǔ)。

3) 推導(dǎo)了高階分?jǐn)?shù)導(dǎo)數(shù)粘彈性模型的數(shù)值求解方法，并應(yīng)用于橡膠隔振元件動力學(xué)模型的數(shù)值求解，通過與試驗結(jié)果進(jìn)行對比，驗證了該方法的正確性與有效性。