甘華權(quán)，李海艷

（廣東工業(yè)大學(xué)機(jī)電工程學(xué)院，廣州 510006）

0 引言

隨著現(xiàn)代科學(xué)技術(shù)的不斷發(fā)展和生活水平的不斷提高，人們對(duì)汽車的功能需求越來(lái)越趨于多樣化和多元化。除了作為日常代步行駛之外，人們對(duì)汽車舒適性和平穩(wěn)性的要求也越來(lái)越高。由于汽車懸掛系統(tǒng)綜合多種作用力，對(duì)車輛行駛的振動(dòng)影響較大，進(jìn)而汽車振動(dòng)直接影響了汽車的穩(wěn)定性、舒適性以及安全性，所以研究汽車懸掛系統(tǒng)的振動(dòng)變形對(duì)改進(jìn)汽車的穩(wěn)定性、舒適性和安全性有著重大的意義[1]。

汽車懸掛系統(tǒng)屬于柔性多體系統(tǒng)的范疇[2]，其在運(yùn)動(dòng)時(shí)會(huì)產(chǎn)生剛性與彈性變形的耦合運(yùn)動(dòng)[3]。所以汽車懸掛系統(tǒng)動(dòng)力學(xué)建模需要考慮彈性變形，但是彈性變形有無(wú)限維自由度的特點(diǎn)，所以無(wú)法求得其精確解。目前比較常用的是通過(guò)有限單元法進(jìn)行離散化處理進(jìn)而求取其近似解[4]。例如，Patil等[5]利用有限元方法來(lái)求解柔性臂運(yùn)動(dòng)到特定位置情況下的撓度。許志華[6]使用有限元方法建立自卸車懸掛系統(tǒng)非線性有限元接觸模型。Nabawy等[7]使用平面桿單元建立了一種綜合的有限元模型來(lái)研究雙叉臂懸架系統(tǒng)的動(dòng)態(tài)響應(yīng)。然而為了方便模型求解，很多學(xué)者使用有限單元法對(duì)汽車懸掛系統(tǒng)動(dòng)力學(xué)建模中一般會(huì)選擇簡(jiǎn)化懸掛系統(tǒng)模型同時(shí)忽略一部分小變形，雖然減少了一些計(jì)算量，但也降低了一些動(dòng)力學(xué)模型的精度，而且有限單元法的精度和網(wǎng)格單元?jiǎng)澐值拇笮∠嚓P(guān)，對(duì)于像汽車懸掛系統(tǒng)這樣的復(fù)雜系統(tǒng)，即使在簡(jiǎn)化了模型的情況下，為了保證計(jì)算精度的網(wǎng)格單元會(huì)劃分的比較小，計(jì)算量也仍然比較大。

因此，本文以汽車懸掛系統(tǒng)為研究對(duì)象，提出了一種利用模態(tài)方法進(jìn)行降維求解的方法，該方法首先通過(guò)拉格朗日方程法建立準(zhǔn)確的汽車懸掛系統(tǒng)的非線性柔性多體動(dòng)力學(xué)模型。然后對(duì)懸掛系統(tǒng)的彈性變形進(jìn)行離散化處理，并通過(guò)有限階級(jí)數(shù)展開(kāi)式精確表示該彈性變形。最后使用廣義α方法迭代求解該模型，并與有限單元法的柔性多體動(dòng)力學(xué)建模求解結(jié)果進(jìn)行對(duì)比分析，得出該方法應(yīng)用于柔性多體動(dòng)力學(xué)建模的優(yōu)勢(shì)。

1 汽車懸掛系統(tǒng)動(dòng)力學(xué)建模

1.1 建立參考坐標(biāo)系

如圖1所示，1、3、5、7為車輪，2、4、6、8為柔性臂，9為車架，10～13為減震器。分析整車懸掛系統(tǒng)需要建立一系列的坐標(biāo)系來(lái)表示其位置，其中oxyz為整體坐標(biāo)系，然后以質(zhì)量塊i(i=1,2,…,9)的重心oi為局部坐標(biāo)系原點(diǎn)，建立與整體坐標(biāo)系平行的初始狀態(tài)的局部坐標(biāo)系oixiyizi。

圖1 汽車懸掛系統(tǒng)

整車懸掛系統(tǒng)根據(jù)解耦控制方法，可以由4個(gè)1/4懸掛系統(tǒng)和一個(gè)為車架與受到的外力組成的系統(tǒng)組成。在整車懸掛系統(tǒng)中，所有的運(yùn)動(dòng)部件全部考慮為均質(zhì)質(zhì)量塊[8]。

如圖2所示，以1/4懸掛系統(tǒng)為例進(jìn)行分析。首先，1為車輪，2為柔性臂，10為減震器，可替換為彈簧阻尼器，質(zhì)量忽略不計(jì)。整體坐標(biāo)系為oxyz，車輪1的局部坐標(biāo)系為o1x1y1z1，柔性臂2的局部坐標(biāo)系為o2x2y2z2（圖中未顯示的x、x1與x2軸為由紙面向里）。由于柔性臂2的長(zhǎng)度遠(yuǎn)大于其截面的直徑，柔性臂的截面方向的變形變化率非常小，所以將柔性臂簡(jiǎn)化為空間梁結(jié)構(gòu)。

圖2 1/4懸掛系統(tǒng)示意圖

假設(shè)在簡(jiǎn)化的空間梁2上存在某一點(diǎn)p，則該點(diǎn)p在局部坐標(biāo)系o2x2y2z2上的位移矢量為：

式中：R2為局部坐標(biāo)系o1x1y1z1的原點(diǎn)在整體坐標(biāo)系中的位置；A2為將柔性臂局部坐標(biāo)系位移轉(zhuǎn)換到整體坐標(biāo)系的旋轉(zhuǎn)矩陣；u2r和u2e分別為p點(diǎn)在局部坐標(biāo)系o2x2y2z2上的剛性和彈性變形位移。

1.2 動(dòng)力學(xué)模型建立

考慮到汽車懸掛系統(tǒng)動(dòng)力學(xué)方程非常復(fù)雜，而使用拉格朗日方程法是直接建立外力與運(yùn)動(dòng)部件變量的關(guān)系，不涉及運(yùn)動(dòng)部件之間的各種約束力，原理相對(duì)簡(jiǎn)單，所以本文利用拉格朗日方程法建立汽車懸掛系統(tǒng)的動(dòng)力學(xué)模型。下面以1/4汽車懸掛系統(tǒng)為例，其拉格朗日方程為[9]：

式中：Qi為運(yùn)動(dòng)部件i廣義坐標(biāo)的廣義外力；λ為拉格朗日乘子向量；Cqi為運(yùn)動(dòng)部件i的約束雅克比矩陣；qi為運(yùn)動(dòng)部件i的廣義坐標(biāo)；Ti為運(yùn)動(dòng)部件i的動(dòng)能。

Ti計(jì)算如下：

式中：q?i為qi對(duì)時(shí)間t的一階微分；Mi為運(yùn)動(dòng)部件i的質(zhì)量矩陣。

根據(jù)式（3）、虛功原理以及邊界條件可求得1/4懸掛系統(tǒng)的動(dòng)力學(xué)方程為：

式中：q?ii為qi對(duì)時(shí)間t的二階微分；Ki為運(yùn)動(dòng)部件i的剛度矩陣；Qvi為運(yùn)動(dòng)部件i二次速度矢量；Qci為運(yùn)動(dòng)部件i的約束雅克比矩陣余量。

其中，可根據(jù)以下式計(jì)算運(yùn)動(dòng)部件i二次速度矢量[10]：

進(jìn)一步的，可以將式（4）寫成以下的矩陣形式：

最后同理，可得整車懸掛系統(tǒng)動(dòng)力學(xué)模型為：

2 動(dòng)力學(xué)模型求解

2.1 動(dòng)力學(xué)方程降維

由于拉格朗日方程法推導(dǎo)得到的汽車懸掛系統(tǒng)柔性多體動(dòng)力學(xué)模型的方程是一組時(shí)變、非線性和強(qiáng)耦合的偏微分-積分方程，所以求汽車懸掛系統(tǒng)動(dòng)力學(xué)模型的解析解非常困難。在這種情況下，一般是將無(wú)限維的動(dòng)力學(xué)模型通過(guò)離散的方法進(jìn)行降維，得到有限維的動(dòng)力學(xué)模型。目前柔性多體動(dòng)力學(xué)模型的常用的降維處理方法主要是有限單元法，該方法通過(guò)劃分單元的大小保證其求解精度。單元?jiǎng)澐衷郊?xì)時(shí)，離散后的動(dòng)力學(xué)模型維度就越高，計(jì)算量就越大，計(jì)算效率就越低。而本文通過(guò)將彈性變形離散成基函數(shù)與對(duì)應(yīng)時(shí)間系數(shù)乘積的級(jí)數(shù)展開(kāi)式進(jìn)而進(jìn)行降維處理，將其代入前面求得的動(dòng)力學(xué)模型即可得到降維后的低維近似模型。

首先根據(jù)柔性臂j的局部坐標(biāo)系ojxjyjzj和柔性臂j的振型函數(shù)取特殊組合的三角函數(shù)為基函數(shù)：

式中：yj為柔性臂j在y軸上的位置表示；Lj為柔性臂j的長(zhǎng)度，(i=1,2,…,N)，ωi為柔性臂的固有頻率。

然后，將柔性臂j的彈性變形離散成式（8）的基函數(shù)與對(duì)應(yīng)時(shí)間系數(shù)乘積的有限級(jí)數(shù)展開(kāi)式表示：

式中：axi(t)、ayi(t)、aji(t)為待求的與時(shí)間相關(guān)的權(quán)函數(shù)。

根據(jù)柔性臂兩端變形為零的邊界條件，可得：

進(jìn)而計(jì)算求得參數(shù)dji為：

最后，本文取前2階級(jí)數(shù)展開(kāi)式精確表示其柔性臂j的彈性變形：

將上述離散的彈性變形代入前面求得的汽車懸掛系統(tǒng)柔性多體動(dòng)力學(xué)模型中，得到降維后的低維柔性多體動(dòng)力學(xué)模型。

根據(jù)上面彈性變形的表示形式可以得到以下剛度矩陣[11]：

由上式可求得柔性臂j的形函數(shù)Sj、幾何矩陣Bj和剛度矩陣Kj，式中qjf為柔性臂j變形對(duì)應(yīng)的廣義坐標(biāo)，E為彈性模量。則有：

式中：Smj為質(zhì)量塊j的截面面積；Igxj為質(zhì)量塊j關(guān)于x方向的慣性矩；Igzj為質(zhì)量塊j關(guān)于z方向的慣性矩。

2.2 動(dòng)力學(xué)方程求解

如圖3所示，通過(guò)拉格朗日方程法可建立的汽車懸掛系統(tǒng)柔性多體動(dòng)力學(xué)模型，通過(guò)求逆矩陣的方法求得初始加速度矢量和拉格朗日乘子。然后將初始時(shí)刻t0時(shí)刻的位移、速度、加速度等矢量作為廣義α方法[12]的初始迭代值，通過(guò)廣義α方法迭代求取t1時(shí)刻的相應(yīng)位移、速度、加速度等，然后以t1時(shí)刻的相應(yīng)位移、速度、加速度為作為t2時(shí)刻廣義α方法迭代的初始值求t2時(shí)刻的相應(yīng)位移、速度、加速度等矢量。然后依次類推，依次求解，直到誤差小于預(yù)設(shè)誤差值且執(zhí)行到未尾時(shí)刻tn則程序迭代結(jié)束，并輸出最后tn時(shí)刻相應(yīng)的位移、速度、加速度矢量等。

圖3 動(dòng)力學(xué)模型求解流程

3 仿真與結(jié)果分析

通過(guò)具體參數(shù)進(jìn)行數(shù)值仿真，車前端、后端的柔性臂長(zhǎng)度分別為0.38 m、0.35 m，車輪、車架質(zhì)量分別為20 kg、300 kg，彈性模量為200 GPa，密度為7800 kg/m3，彈簧彈性系數(shù)為10000 N/m，車架質(zhì)心( 0.25,0,0.36)m，柔性臂截面內(nèi)、外圓半徑分別為0.015 m，0.01 m。

圖4所示為沿x軸方向的彈性位移的比較，實(shí)線為本文方法求解得到的x軸彈性位移，虛線為有限元方法求解得到的x軸彈性模型，第一到第四行分別為柔性臂2、4、6、8的x軸彈性位移對(duì)比。

圖4 不同時(shí)刻柔性臂x軸變形位移

圖5所示為柔性臂中點(diǎn)沿z軸方向的彈性位移的比較，其中實(shí)線為本文方法表示的動(dòng)力學(xué)模型求解得到的z軸彈性位移，虛線為有限元方法表示的動(dòng)力學(xué)模型求解得到的z軸彈性模型，圖5（a）～（d）分別為柔性臂2、4、6、8的z軸彈性位移對(duì)比。

圖5 不同柔性臂中間節(jié)點(diǎn)變形位移

從圖4、圖5可知：首先，將彈性變形降維到64維的動(dòng)力學(xué)建模求解的有限元方法與將彈性變形降維為二維的動(dòng)力學(xué)模型求解的本文方法計(jì)算得到的不同時(shí)刻柔性臂x軸變形位移及不同柔性臂中間節(jié)點(diǎn)變形位移趨勢(shì)一致，且結(jié)果相差無(wú)幾，即精度相當(dāng)；其次，本文方法求解所用時(shí)間僅為0.79 s，而有限元方法求解所用時(shí)間為2.3 s，其所用時(shí)間是本文方法的4倍多，即本文方法的計(jì)算效率遠(yuǎn)高于有限單元方法；最后，如果進(jìn)一步提升求解動(dòng)力學(xué)模型的精度，有限單元法的網(wǎng)格劃分將會(huì)更多，計(jì)算量將大大增大，而本文則是只需增加一兩階基函數(shù)即可，所以在相同精度下，本文方法在求解效率方面優(yōu)勢(shì)非常顯著。

4 結(jié)束語(yǔ)

本文通過(guò)以整車汽車懸掛系統(tǒng)作為研究對(duì)象，考慮彈性變形對(duì)汽車懸掛系統(tǒng)的影響，將柔性臂簡(jiǎn)化成柔性空間梁結(jié)構(gòu)，不忽略柔性臂的軸向和截面變形，使建立的汽車懸掛系統(tǒng)的柔性多體動(dòng)力學(xué)模型比以往的將柔性臂簡(jiǎn)化成平面梁表示動(dòng)力學(xué)模型更完整。通過(guò)對(duì)比整車汽車懸掛系統(tǒng)動(dòng)力學(xué)模型有限單元法降維求解和本文所提方法可得，兩種方法都可以保證整車汽車懸掛系統(tǒng)動(dòng)力學(xué)建模的精度，在求解精度相差不大的情況下，本文所提方法的計(jì)算效率明顯要高得多。

為了進(jìn)一步研究汽車懸掛系統(tǒng)運(yùn)動(dòng)中產(chǎn)生的彈性變形對(duì)汽車懸掛系統(tǒng)振動(dòng)的影響，以后的研究通過(guò)建立更高維的柔性多體動(dòng)力學(xué)模型和研究其振動(dòng)頻率進(jìn)行分析。