基于譜方法的柔性懸架系統高效動力學仿真*

李曉飛,李海艷,梁桂銘

(廣東工業大學 機電工程學院,廣東 廣州 510006)

0 引 言

車輛懸架系統受到多種作用力,其對車輛系統動力學性能影響非常大,所以建立懸架系統動力學模型并進行仿真是研究車輛系統動力學性能的重要工作之一。

車輛的懸架系統屬于柔性多體系統范疇,其柔性梁在運動時存在大位移的剛性運動,同時產生一定的彈性變形運動,這兩種運動是高度耦合的[1]。在懸架系統的動力學建模中,由于彈性變形有無限維自由度的特點,其精確解無法求得。目前處理彈性變形的方法一般是將其在時間和空間方面離散,求取其近似解[2]。常用的離散方法有3種,即有限單元法、假設模態法和譜方法。

其中,有限單元法是在空間上將柔性體劃分成網格單元,通過求解網格單元的變形近似彈性變形的精確解。PATIL A M等人[3]使用有限元方法,建立了懸架系統靜止狀態下的動力學模型。GADADE M B等人[4]針對商用車前下懸架臂的設計問題,研究了靜止狀態下的懸架系統三維有限元動力學模型。NABAWY A E等人[5]使用平面桿單元,建立了一種懸架系統綜合的有限元模型,研究了雙叉梁懸架系統的動態響應。

依據有限單元法離散求解了柔性體彈性變形,雖然保證了求解精度,但是其劃分網格單元較多,使得其計算量非常大,計算效率低。

而采用假設模態法時,需通過先驗知識求得柔性體的模態和振型,然后根據精度要求決定截斷振型函數的階數,再通過振型擬合柔性體的變形。VAKIL M等人[6]使用假設模態法建立了柔性連桿動力學模型,并分析了其動力學模型的可行性。潘云[7]采用假設模態法,建立了柔性機械臂的動力學模型,并對機械臂的軌跡控制進行了研究。張曉宇等人[8，9]依據假設模態法,建立了柔性連桿的動力學模型,研究了其對控制器設計的影響。

雖然采用假設模態法擬合彈性變形函數比較簡便,但是其所增加的振型函數的項數會產生計算效率急劇降低的問題。

譜方法[10]是通過一組正交多項式的有限級數來表示柔性體彈性變形。它無需對網格進行劃分,也不需要先驗知識處理彈性變形函數的擬合,具有收斂快、高精度的特點。

基于上述原因,筆者將譜方法引入到柔性多體動力學建模中,將巴哈賽車懸架系統中柔性梁的彈性變形函數離散成Chebyshev多項式的有限級數形式,通過置換方法處理彈性變形函數,再利用彈性變形函數和拉格朗日動力學公式建立懸架系統的柔性多體動力學模型,通過廣義α方法[11]進行迭代求解;最后將計算結果與其他方法進行對比,驗證基于譜方法的懸架系統動力學建模求解的準確性和高效性。

1 懸掛系統解耦與柔性梁運動描述

巴哈賽車模型如圖1所示。

圖1 巴哈賽車模型1,3,5,7—輪胎;2,4,6,8—柔性梁;9—車架;10,11,12,13—減震器

為了便于對懸架系統進行建模,筆者依據浮動框架建立懸架系統的坐標系,即以初始狀態車架的重心o為原點建立oxyz絕對坐標系,然后以運動部件i(i=1,2,…,9)的重心oi為局部坐標系原點,建立初始狀態與絕對坐標系平行的局部坐標系oixiyizi;

由于懸架系統左右對稱,筆者依據文獻[12]的懸架解耦控制方法,將懸架系統解耦成4個1/4懸架和1個受到4個彈簧阻尼力的車架系統,整個懸架系統按5部分進行動力學建模,再利用其作用力之間的關系形成完整的動力學模型。

接下來,筆者對1/4懸架系統進行分析。

1/4懸架系統原理圖如圖2所示。

圖2 1/4懸掛系統k1,k10—彈簧阻尼器1,10彈性系數;c1,c10—彈簧阻尼器1,10彈性系數

筆者將減震器簡化成彈簧阻尼器10,將輪胎1簡化成剛體1和彈簧阻尼器1。由于柔性梁2的長度遠大于其截面的直徑,柔性梁的截面方向的變形變化非常小,其中心線的變形可比較精確地表示截面的變形,柔性梁可簡化為空間梁結構。

筆者依據圖2對柔性梁2進行運動分析。

柔性梁2上某一點p的位移矢量為:

r2=R2+A2u2=R2+A2(u20+u2f)

(1)

式中:R2—點o2在絕對坐標系oxyz上的位移;A2—柔性梁2的旋轉矩陣;u20,u2f—點p在局部坐標系o2x2y2z2上的剛性位移和彈性變形位移。

2 基于譜方法的動力學建模

2.1 柔性梁彈性變形表示

根據柔性梁2的運動描述可知,柔性梁2的位移由剛性位移與彈性變形位移疊加而成。因此,柔性梁2的彈性變形位移可表示為:

(2)

式中:u2fx,u2fy,u2fz—彈性變形在x2,y2和z2軸上的投影。

筆者依據譜方法將彈性變形位移函數進行離散,將空間項和時間項分離,得到有限級數表達式:

(3)

式中:Tm(y)—區間[-L2/2,L2/2]上的第一類切比雪夫多項式;a2m(t),b2m(t),c2m(t)—對應的譜表示系數;L2—柔性梁2長度;M—多項式截斷數。

將彈性變形函數寫成矩陣表示形式,即:

u2f=U2fqU2+V2fqV2

(4)

其中:

(5)

(6)

(7)

由于根據譜方法得到的柔性梁2彈性變形位移函數需要滿足相應的邊界條件,需要對彈性變形位移函數進行置換處理。

柔性梁2上點A為固定連接,點B為鉸接,所以u2fx,u2fy,u2fz需要滿足如下邊界條件:

(5)

將式(8)寫成矩陣形式,即:

U2qU2+V2qV2=0

(9)

其中:

(10)

根據式(9),qV2可以表示為:

(11)

將式(11)中的qV2置換代入式(4),可得:

(12)

式中:u2f—滿足邊界條件的譜表示彈性位移函數;S2f—譜表示彈性位移函數的形函數;qU2—譜表示彈性位移函數的廣義坐標。

2.2 動力學模型

根據式(1,12)所求得的柔性梁2的譜表示彈性位移函數u2f,形函數S2f和彈性變形廣義坐標qU2,可得柔性梁2的譜表示的位移矢量為:

r2=R2+A2(u20+S2fqU2)

(13)

依據位移矢量r2,對t求微分得到速度矢量,可將柔性梁2的動能[13]表示為:

(14)

式中:ρ2—柔性梁2的密度;M2—質量矩陣;V2—柔性梁2的三維可行域;“·”—對時間t微分。

其中:

(15)

式中:θ2—局部坐標系o2x2y2z2中x2、y2、z2軸與絕對坐標系oxyz中x、y、z軸的夾角向量。

將式(14)代入拉格朗日方程可得:

(16)

經整理可得譜表示的柔性梁2的運動方程為:

(17)

式中:Cq2—雅克比矩陣。

其中:

(18)

(19)

式中:D—三維彈性矩陣。

因為柔性梁2所受外力只有重力,所以廣義外力Q2e可通過虛功原理[14]求得,即:

(20)

式中:g—重力加速度;δq2—虛位移。

由于車輪1簡化為剛體1與彈簧阻尼器1,其剛體1運動方程為:

(21)

式中:M1—剛體1的質量矩陣;R1—剛體1的質心坐標。

其中:

(22)

式中:θ1—局部坐標系o1x1y1z1中x1、y1、z1軸與絕對坐標系oxyz中x、y、z軸的夾角向量。

由于剛體1所受的外力為彈簧阻尼力和重力,且作用于質心,其廣義外力Q1e為:

(23)

其中:

(24)

式中:γ—彈簧阻尼器10與z軸的夾角;u9z—車架的在z9軸的位移;β9—車架關于y9軸的轉角;La—車架質心在x軸向與彈簧阻尼器10的距離;Rf—地面激勵。

根據式(17,21)可得1/4懸架系統動力學方程為:

(25)

由于廣義坐標q1與q2受約束影響而不完全獨立,Cq12需要根據柔性梁2與車輪1的約束方程求得。

柔性梁2與車輪1在點A為固定連接,在點B與底盤鉸接,其對應的約束方程為:

(26)

式中:u1A,u20A—點A在局部坐標系o1x1y1z1和o2x2y2z2上的初始位移;S2fA—點A在局部坐標系o2x2y2z2上的形函數;S2fB—點B在局部坐標系o2x2y2z2上的形函數;RcB—點B在全局坐標系oxyz上的位移。

將式(26)表示為:

Φq12=0

(27)

雅克比矩陣Cq12和加速度右項Qc12可由文獻[15]中的公式求得,即:

(28)

式中:“‥”—對t的二階微分。

由于車架與車輪是通過彈簧阻尼器連接,車架運動只與外力有關,且車架為剛體,其動力學方程為:

M9q9=Q9e

(29)

式中:M9—車架的質量矩陣;q9—廣義坐標。

由于懸架系統[16]左右對稱,車架受到的外力為彈簧阻尼力和重力,Q9e為:

(30)

其中:

(31)

(32)

式中:Lb—車架質心在x軸向與彈簧阻尼器11的距離;k11,k12,k13—彈簧阻尼器11-13的彈性系數;c11,c12,c13—彈簧阻尼器11-13的阻尼系數;ujz—車輪j(j=1,3,5,7)的z軸的位移。

同理,通過式(25),可求其他1/4懸架系統的動力學方程,然后根據所求得的懸架系統動力學方程和車架的動力學方程,依據廣義坐標,可形成整車懸架系統的動力學模型。

將其與約束方程整理成矩陣形式為:

(33)

3 數值仿真及分析

3.1 廣義α方法迭代求解

懸架系統參數[17]如表1所示。

表1 懸架系統參數

為了驗證基于譜方法的懸架系統動力學模型的可行性和高效性,筆者利用表1中的參數對其進行仿真求解。

通過廣義α方法的仿真求解流程如圖3所示。

圖3 仿真流程圖

具體仿真流程為:

由于流程圖3中其他參數和變量相應文獻有詳細描述和具體計算公式,此處不再贅述。

3.2 與有限元動力學模型的車架位移對比

筆者將仿真的車架位移與求解三維有限元動力學模型得到的車架位移相對比。

其中,車架z軸位移如圖4所示。

圖4 車架z軸位移

車架y向轉角如圖5所示。

圖5 車架y向轉角

通過二范數可計算其相對誤差ζ:

(34)

式中:DF—精確解;DS—試驗解。

以有限元的仿真的車架位移為精確解DF,基于譜方法動力學模型的車架位移為試驗解DS,車架z向位移二范數相對誤差為3.57%,y向轉角的二范數相對誤差為3.37%。

該結果表明:基于譜方法的懸架系統動力學模型的誤差在可容許范圍內。

3.3 與有限元動力學模型柔性梁彈性位移對比

根據仿真迭代求出的廣義坐標q和已知的彈性變形的形函數Sjf(j=2,4,6,8),可求出柔性梁j的彈性變形。

筆者將不同Chebyshev多項式階數下,基于譜方法的動力學模型仿真求得的柔性梁彈性變形位移,與有限元動力學模型求解得到的不同時刻下彈性變形位移相對比,結果如圖(6～11)所示。

圖6 柔性梁2的x向彈性位移

圖7 柔性梁2的y向彈性位移

圖8 柔性梁2的z向彈性位移

圖9 柔性梁4的x向彈性位移

圖10 柔性梁4的y向彈性位移

圖11 柔性梁4的z向彈性位移

接著,筆者以有限元動力學模型求解出的柔性梁彈性位移為準確值,以基于譜方法動力學模型計算出的柔性梁彈性變形為試驗值,通過式(34)計算彈性變形二范數相對誤差,結果如表2所示。

表2 不同階數下的柔性梁相對誤差

根據表2的數據和圖(6～11)的結果可得:

(1)隨著譜方法選取Chebyshev多項式項數的增加,筆者所提的方法的計算精度不斷提高,當所取Chebyshev多項式項數為24時,柔性梁彈性變形最大相對誤差為1.18%;

(2)隨著所取Chebyshev多項式項數增加,相對誤差最大值的遞減幅度越來越小,當所取多項式項數達到16時,繼續增加多項式項數相對誤差降低幅度已經不顯著,最大降低幅度為0.54%;

(3)表2所示的時間為在相同計算機軟件與硬件條件下,懸架系統動力學模型的求解時間;依據表2中的數據可知,隨著所取Chebyshev多項式項數遞增,求解總時間也在依次遞增;

(4)當振型函數和Chebyshev多項式都取24階時,假設模態法動力學模型求解時間是基于譜方法的懸架系統動力學模型求解計算時間的4.896倍。

4 結束語

筆者采用譜方法將懸架系統中柔性梁的彈性變形函數離散成Chebyshev多項式的有限級數形式,通過置換方法處理彈性變形函數得到的邊界條件,再利用彈性變形函數和拉格朗日動力學公式建立懸架系統的柔性多體動力學模型,通過廣義α方法迭代求解,最后將基于譜方法的懸架系統動力學模型的求解結果與其他方法下懸架系統動力學模型求解結果進行對比,得到了以下結論:

(1)將基于譜方法的懸架系統動力學模型的求解結果與有限元動力學模型仿真的求解結果相對比,車架位移的相對誤差為3.57%,柔性梁的彈性變形相對誤差最小為1.18%;

(2)當振型函數和Chebyshev多項式都取24階時,基于譜方法的懸架系統動力學模型的計算效率是假設模態法動力學模型計算效率的4.896倍,體現了譜方法的高效性;

(3)當所取Chebyshev多項式項數達到16時,繼續增加階數,基于譜方法的懸架系統動力學模型的計算精度提升不大,計算效率降低比較多,所以在動力學建模中需要適當地選取Chebyshev多項式項數。

目前譜方法使用的正交多項式是考慮截斷其前M項,不是有選擇性地挑選正交多項式中的項。因此,在后續的研究中,筆者將通過壓縮感知稀疏理論,有選擇性地挑選正交多項式中項,來擬合待求函數,以進一步提升其計算精度和效率。