農用車輛橡膠隔振元件動力學模型研究*

高琦，邱緒云，宋裕民

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

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

0 引言

我國農業機械化迅猛發展的過程中，農用車輛的廣泛應用大大提高了農業生產和運輸的效率[1]。農用車輛的使用工況相對惡劣，在作業和行駛過程中使車輛各結構部件產生劇烈振動，同時發動機產生的振動也會傳遞給車輛各系統及其部件，不但會降低農用車輛作業穩定性，引起結構噪聲，也會影響各部件的疲勞壽命[2-3]。

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

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

1 橡膠隔振元件動態力學試驗

選取某農用車輛發動機前懸置軟墊為試驗對象，基于MTS831襯套試驗機進行了軸向動態加載試驗，如圖1所示。

圖1 某農用車輛發動機前懸置軟墊軸向動態加載試驗Fig. 1 Axial dynamic loading test for the engine front mounting cushion of an agricultural vehicle

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

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

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

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

2 橡膠隔振元件動力學模型

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

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

2.1 彈性單元

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

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

(1)

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

x——加載位移，mm；

a0,a1,a2,a3,…,an——冪次項的系數。

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

(2)

2.2 摩擦單元

摩擦單元中，摩擦力與加載位移的關系[13]

(3)

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

Ffmax——最大摩擦力，N；

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

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

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

(xs，Ffs)——摩擦力與位移變化曲線上的狀態參考點，此處取(0，0)。

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

(4)

一個加載循環損耗的能量為

(5)

其中u0=Ff0/Ffmax。

2.3 粘彈性單元

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

圖5 高階分數導數粘彈性模型示意圖Fig. 5 High-order fractional derivative viscoelastic model

根據圖5所示模型，可得粘彈力與位移的關系

λ1λ2x(t)

(6)

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

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

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

α、β、γ——取值范圍(0，1)的分數導數階數，其中α+γ的最高階數為2。

對式(6)進行傅里葉變換，得該分數導數粘彈力模型的復剛度

(7)

(8)

(9)

一個加載循環損耗的能量為

(10)

2.4 橡膠隔振元件動力學模型

將彈性單元、摩擦單元與分數導數粘彈性單元疊加，得到橡膠隔振元件動力學模型的總響應力與加載位移的關系

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

(11)

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

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

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

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

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

(12)

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

(13)

橡膠隔振元件模型的阻尼系數

(14)

3 橡膠隔振元件動力學模型參數辨識

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

3.1 摩擦單元參數辨識

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

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

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

(15)

3.2 彈性單元參數辨識

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

3.3 粘彈性單元參數的辨識

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

粘彈性單元參數辨識實質上是一個最優參數估計問題，為使識別結果盡可能接近試驗測量結果，通過建立橡膠元件模型計算的動態響應與試驗得到元件的動態響應誤差最小的優化函數來實現。此處基于橡膠隔振元件動剛度和阻尼系數誤差建立優化目標函數

(16)

式中：M——辨識工況的總個數；

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

以式(16)所示目標函數作為適應度函數，算法中更新粒子的方法如式(17)所示。

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

(17)

式中：t——當前迭代次數；

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

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

bi(t)——一個隨機產生的概率；

(18)

式中：d——粒子的搜索維數；

n——群中的粒子數；

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

(19)

式中：ωEt——動態調整的慣性權重；

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

κ——社會影響系數，κ取值比較大時，可能導致過早收斂到平均行為的現象；

r1、r2——廣義多維學習部分和社會學習部分的權重，取值為0～1之間均勻分布的隨機數。

(20)

式中：ni——落入每個小區間的粒子數。

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

(21)

參數識別過程中算法的收斂性準則為

(22)

式中：σfobj——所有粒子適應值的標準差；

μfobj——所有粒子適應值的平均值。

該判據可以避免在fobj的全局最優值附近出現小的波動。否則，當迭代次數達到最大迭代次數Gmax時，優化過程終止。

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

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

3.4 橡膠隔振元件模型參數辨識結果與驗證

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

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

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

表1 橡膠隔振元件模型各單元參數辨識結果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 橡膠隔振元件阻尼系數模型計算值與試驗值對比Fig. 9 Comparison of damping coefficients between model calculation value and test value of rubber vibration isolation element

4 橡膠隔振元件動力學模型數值求解

為了將建立的橡膠隔振元件動力學模型應用于農用車輛動力學建模與分析中，對橡膠隔振元件模型的數值求解方法進行推導。

4.1 高階分數導數粘彈性單元數值求解

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

(23)

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

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

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

(24)

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

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

(25)

由于

(26)

因此，可得到

(27)

零初始條件下，當N取一極大值時，可進一步得到分數階導數的數值計算公式

(28)

(29)

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

考慮如下形式的分數階微分方程

0

(30)

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

Fs(t)——響應力，N；

sj——分數階導數項的系數；

sl——一次項的系數；

σl、σl-1——分數階導數的階數，σl和σl-1為大于0的實數，且σl>σl-1。

式中：θh——取值在(0，1]區間內的實數。

根據分數導數運算的可加性，對式(30)進行降階處理，設

(31)

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

(32)

(33)

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

4.2 數值求解結果驗證

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

分別取頻率1 Hz、振幅2 mm和頻率50 Hz、振幅1 mm兩種正弦信號作為位移輸入，進行數值求解得到橡膠隔振元件的總響應力，并與試驗結果進行對比，如圖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可以看出，在不同頻率工況下，數值求解得到橡膠隔振元件的力—位移曲線與試驗測得的力—位移曲線具有較好的一致性。

5 結論

1) 以某農用車輛發動機前懸置軟墊為載體，通過動態加載試驗對其軸向頻率響應特性和幅值響應特性進行了研究，建立了彈性單元、摩擦單元和高階分數導數粘彈性單元疊加的橡膠隔振元件動力學模型。

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

3) 推導了高階分數導數粘彈性模型的數值求解方法，并應用于橡膠隔振元件動力學模型的數值求解，通過與試驗結果進行對比，驗證了該方法的正確性與有效性。