剛性輪—土壤相互作用的多體動力學仿真研究

摘 要:車輪與松軟地面相互作用的研究對越野車輛的設計、性能評價和行駛操作等至關重要。提出了一種新的多體系統大規模接觸碰撞理論與地面力學理論相結合的車輪土壤相互作用仿真方法。基于自行開發的多體動力學仿真軟件研制了車輪土壤相互作用模塊并構建了單輪土槽多體仿真模型。對相同實驗條件下得到的仿真與實驗結果的對比分析表明，該仿真算法具有較高的工程實用性。

關鍵詞:車輪土壤相互作用; 多體動力學; 接觸碰撞; 車輛地面力學; 單輪土槽

中圖分類號:TP319.9文獻標志碼:A

文章編號:1001-3695(2009)09-3289-04

doi:10.3969/j.issn.1001-3695.2009.09.025

Dynamic multibody simulation research of rigid wheel-soil interaction

ZUO Yan-rui ZONG Zhi-jian LIU Zhong-tu CHENG Yuan FAN Shi-chao2

(1. Advanced Design Technology Center， College of Engineering， Sun Yat-sen University， Guangzhou 510006， China; 2. Beijing Institute of Spacecraft Environment Engineering， Beijing 100094， China)

Abstract:The study of wheel-soft terrain interaction is important for design， performance evaluation and steering control of off-road vehicle. This paper proposed a new method to realize the wheel-soil interaction simulation which applied the large-scale multibody contact and terramechanics theory. Based on self-developed multibody dynamics simulation software， it deve-loped a simulation module and constructed a multibody model of the single wheel-soil testbed. The comparison between simulation and experimental results obtained under the same experimental conditions proves the method’s practical applicability.

Key words:wheel-soil interaction; multibody dynamics; contact problem; vehicle terramechanics; single wheel-soil testbed

0 引言

越野行駛中，松軟地面作用于車輪的滾動阻力、地面推力和有效轉動力矩等是影響越野車輛行駛性能的直接因素，開展車輪土壤相互作用課題的研究對越野車輛的設計、性能分析、整車作業效率和松軟地面通過性評價等至關重要。至今為止，這方面的研究工作已有不少[1~3]。特別是20世紀70年代以來，隨著計算機技術和計算技術的發展，在車輪土壤相互作用課題研究中涌現了大量采用計算機仿真技術得到的研究成果。例如，Andrade等人[4]用有限平面元素方法建立了沙土模型;Wanjii等人[5]使用土壤的粘彈性模型建立了車輪下面土壤的應力分布模型，預測了車輪的行駛性能;Reina等人[6]在此模型基礎上提出了改進并測量了車輪的滑轉率和沉陷量;Nakashima等人[7]用有限元與離散元相結合的方法對車輪和土壤的相互作用進行了研究。然而，這些研究成果由于使用不便且通常只能計算某種力或某種效果，以及經常需要很長的計算時間(如有限元法)等原因，難以滿足越野行駛仿真的工程需要。

注意到車輪土壤相互作用的本質是車輪與大量土壤微元的大規模接觸碰撞問題，以及計算多體系統動力學仿真軟件所具備的優點，即可以最大限度地開發計算機對室內土槽系統、車輛行走機構等復雜系統運動學、動力學和控制系統的分析與綜合能力。筆者提出了一種新的研究車輪土壤相互作用問題的數值仿真方法。該方法通過三角面網格技術與車輛地面力學理論相結合，建立土壤離散多體力學模型，通過車輪與大規模土壤離散單元的接觸碰撞理論和數值求解技術，仿真車輪與土壤間的相互作用。根據此仿真方法與思路，筆者基于自己開發的多體動力學仿真平臺，研制了車輪土壤相互作用仿真模塊，并構建了實驗室單輪土槽系統的多體動力學仿真模型，進行了車輪受力、車輪沉陷、車輪滑轉等重要內容的仿真研究。

本文的討論重點是車輪土壤相互作用力的仿真計算。首先建立了土壤離散多體模型和車輪土壤系統模型，給出了仿真車輪土壤相互作用的大規模接觸碰撞動力學方程。通過構建實驗室單輪土槽裝置的多體動力學仿真系統，仿真了車輪不同滑轉率、不同載荷等工況下的運動并輸出了車輪土壤相互作用力仿真結果，最后進行了與實驗室土槽測量實驗結果的誤差對比與分析。

1 土壤多體動力學建模

土壤具有復雜的物理特性和力學特性，車輛地面力學仿真研究工作者提出了多種土壤建模思想和方法[7~9]。根據工程實際需要對仿真速度和仿真效果提出的要求，本文基于文獻[9]的建模思想，提出了采用面三角網格將松軟土體分割成大規模均勻且緊密排列的離散單元體的方法，如圖1所示。其中，單元網格尺寸在[1~10 mm]內取值。每個單元體的面單元節點代表由此單元網格里的所有土壤顆粒組成的土壤顆粒簇。此方法不僅可以減少計算量、縮短計算時間、降低對計算機性能的需求，同時可以兼顧仿真結果的真實性。

鑒于土壤復雜的重復加載—卸載—沉陷特性[10]以及在平板沉陷實驗中呈現出的土壤應力隨平板穿透速度增加而增加的阻尼特性[11]，本文參考文獻[12]，采用圖2所示的粘彈塑性力學模型建立土壤離散多體動力學模型，以充分反映土壤的粘彈塑特性。圖中虛線框部分表示用剛度系數cf和阻尼系數η描述的土壤粘彈性特性，實線框部分表示用屈服常量μp和壓縮系數μv描述的土壤塑性特性。本模型描述的土壤動態承載—沉陷關系如下:

p(t)=cf(z(t)-zp(t))+η(z#8226;(t)-z#8226;p(t))(1)

p(t)=μvzp(t)+μp(2)

對式(1)(2)求導并整理得

(cf+μv)p(t)+ηp#8226;(t)=cfμvz(t)+ημvz#8226;(t)+cfμp(3)

式(1)～(3)中，z為土壤在車輪荷載作用下的總沉陷量，zp為土壤塑性沉陷量。其中，式(1)描述了土壤粘彈性沉陷平衡關系;式(2)描述了土壤塑性沉陷平衡關系。不同于Bekker靜態承壓—沉陷半經驗公式[1]，本模型全面描述了土壤特性，是土壤承壓—沉陷的細觀、動態模型，非常適用于車輪土壤相互作用的實時仿真。

2 車輪土壤系統接觸碰撞模型的建立

在車輪土壤相互作用仿真的初步工作中，不考慮車輪的轉向問題，車輪土壤系統可簡化為如圖3所示的平面系統。在這個平面系統中，只需考慮作用于車輪的平面力和垂直于此平面(xoz平面)的轉動力矩。

車輪運動過程中，仿真引擎通過接觸碰撞檢測模塊追蹤并記錄車輪土壤的接觸碰撞信息。設某tc時刻碰撞檢測模塊檢測到車輪與N個土壤單元發生碰撞，與M個土壤單元發生接觸，則應用第1章中建立的土壤離散多體動力學模型，建立車輪土壤系統的大規模接觸碰撞動力學模型，如圖3所示。其中:V為車輪質心速度;ω車輪轉動角速度;θf為車輪前進角;θr車輪離去角;θi為第i個車輪土壤接觸點與車輪豎直軸線的夾角;W為車輪有效載荷;M為驅動力矩;(Fi)ni為第i個車輪土壤接觸點處的二維平面接觸力的法向分力;(Fi)ti為第i個車輪土壤接觸點處的車輪土壤二維平面接觸力的切向分力。

針對tc時刻檢測到的車輪土壤接觸碰撞信息，仿真引擎的接觸碰撞響應模塊的計算步驟如下:

a)對N個碰撞進行響應，計算出N個碰撞沖量并將其應用于車輪，更新車輪碰撞后的運動狀態。

b)重新檢測碰撞響應前記錄的M個車輪土壤接觸點處的法向相對加速度(ai)ni(i=1，2，NA1AD，M)。如果檢測到某些接觸點處有(ai)ni>0(i=1，2，NA1AD，M)，則認為車輪在這些接觸點處即將與土壤分離，不再需要對其進行接觸力的計算;反之，則記錄接觸點的索引、位置、法線等相關信息。

c)所有M個接觸檢測完畢后，假設車輪與土壤間依然存在m個接觸，則參考文獻[13]建立車輪土壤系統的接觸動力學模型以及線性接觸補償條件:

a=AF+b，withan(tc)≥0Fn(tc)≥0Fn(tc)#8226;an(tc)=0(4)

其中，廣義相對加速度向量a∈R(2m)×1和廣義接觸力向量F∈R(2m)×1的具體形式如下:

a=((a1)n1，(a1)t1，(a2)n2，(a2)t2，NA1AD，(am)nm，(am)tm)TF=((F1)n1，(F1)t1，(F2)n2，(F2)t2，NA1AD，(Fm)nm，(Fm)tm)T

上述為一個線性補償問題，應用求解這類問題的成熟的Baraff數值算法[13]求解得到包含車輪土壤系統所有接觸點處接觸力的廣義接觸力向量F。

由于車輛在松軟地面上行駛時經常發生打滑現象，本文將車輪土壤的接觸問題歸結為動摩擦接觸問題。根據描述松軟土壤弱抗剪特性的Janosi公式[14]，并結合描述物體動摩擦接觸問題中切向接觸力與法向接觸力之間關系的庫侖摩擦定律，本文參考文獻[15]提出了如下適用于本仿真的切向摩擦力計算模型:

(Fi)ti=min(cSi+(Fi)nitan φ，μt(Fi)ni)#8226;(1-e-ji/k)(5)

其中:c為土壤內聚力;φ為土壤內摩擦角;Si為接觸點i所在接觸面元的面積;μt為動摩擦系數;k為土壤水平剪切位移模量;ji為車輪上接觸點i的剪切位移，與其所在位置θi、車輪前進角θf和車輪滑轉率s有關。經推導，有如下計算式:

ji=r(θf-θi)-(1-s)(sin θf-sin θi)」(6)

3 車輪受力計算模型

計算得到車輪土壤接觸區域里的廣義接觸力向量F后，對圖3中所示的接觸區域里所有接觸點處的法向接觸力(Fi)ni的水平分量求和，得到車輪的滾動阻力:

R=∑mi=1[(Fi)ni#8226;sin θi](7)

對所有接觸點處的切向接觸力(Fi)ti的水平分量求和，得到車輪的地面推力:

H=∑mi=1[(Fi)ti#8226;cos θi](8)

得到車輪的滾動阻力和地面推力之后，即有車輪掛鉤牽引力:

DP=H-R(9)

所有切向接觸力(Fi)ti對車輪輪心的力矩之和就是車輪的有效轉動力矩:

Me=∑mi=1[r×(Fi)ti](10)

4 土槽仿真系統

帶有精密測控設備的室內土槽測試系統是地面力學中的一種重要實驗設備，使用土槽測試系統可以測量行走機構的滾動阻力、提供的牽引力、行駛的滑轉率、牽引效率和車輪在土壤中的下陷量等參數。土槽實驗不受天氣等因素的影響，土壤參數和加載條件易于控制，其實驗結果具有可對比性[16]。

為檢驗車輪土壤相互作用仿真模塊的仿真結果，本文在自行研制的多體系統動力學仿真平臺上構建了仿真界面如圖4所示的土槽仿真系統。

4.1 土槽仿真平臺構建

單輪土槽仿真系統的機械部件如圖4所示，主要由三大部分組成:a)控制車輪運動的導軌、臺架和車輪軸;b)車輪;c)土槽。各部件間的物理約束情況如表1所示。部件間的相對運動通過驅動約束控制，如通過圓柱副—距離驅動約束控制臺加速度為0.03 m/s，通過角驅動約束控制車輪轉動角速度在0.1~0.4 rad/s范圍內取值等。

表1 仿真系統物理約束

約束被約束剛體約束允許的剛體運動

固定約束導軌無

圓柱副—距離驅動約束臺架與導軌臺架和導軌可沿圓柱副的z軸(與世界坐標系的x軸方向相同)轉動和移動，其中移動速度可控

圓柱副約束車輪基座與臺架車輪基座與臺架可沿圓柱副的z軸(與世界坐標系的z軸方向相同)轉動和移動

轉動副約束車輪與車輪基座車輪與車輪基座可沿轉動副的z軸(與世界坐標系的y軸方向相同)轉動

角驅動約束車輪與導軌驅動車輪轉動，并控制車輪的轉動角速度

4.2 被試車輪和土壤參數

影響車輪土壤相互作用效果的車輪自身設計參數和土壤力學參數很多，其中剛性車輪的荷載、直徑、胎面寬度、胎面形狀和材料、輪刺高度、輪刺數量等參數和土壤的沉陷指數n、變形模量k、切線模量K、內聚力c、內摩擦角φ、壓縮系數μv、屈服強度μp、土壤剛度系數cf和阻尼系數η等參數是影響車輪受力的重要參數。表2列寫了被試車輪和土壤的信息。

4.3 實驗方案設計

為進行仿真結果的對比分析及驗證，本文依據三種車輪載荷，多個滑轉率正交結合的實驗方案，進行了多組仿真實驗。其中，車輪載荷的選取以質量輕、工作環境重力場小的月球車車輪荷載為依據，分別選取W1=100 N，W2=200 N和W3=300 N。由車輛地面力學實驗經驗，車輪最佳滑轉率大約在10%，因此在滑轉率選取上采用0~10%漸密，10%左右最密，10%~100%漸疏的原則。

5 結果分析

為進行仿真結果的有效分析，筆者在北京航空航天大學汽車工程系設計的行星探測車車輪行走性能實驗臺上同時進行了測量實驗。該實驗臺主要包括土槽體、機械系統、數據采集與控制系統。測量與仿真的實驗條件以及被試車輪、土壤等完全對應。

在土槽仿真系統的仿真運算過程中，取仿真時間步長Δt=0.01 s。對大量仿真實例計算時間的統計顯示，仿真系統30 s的運動，普通配置臺式計算機運行的時間大約為20 min左右。隨著計算機硬件性能的發展以及本仿真軟件計算技術的進一步改進，提高仿真速度，實現實時仿真是可能的。以車輪在W=300 N荷載下，以角速度ω=0.196 rad/s運行時的一個具體仿真實例為例，繪制圖5以說明本仿真的輸出情況。圖中以驅動力矩、滾動阻力兩個指標為例，橫軸為隨意選取的一段仿真時間點。選取車輪運動穩定后的18~23 s的仿真輸出結果，去除個別噪點，繪制車輪驅動力矩和滾動阻力隨時間的變化曲線，如圖6和7所示。兩張圖片均顯示出車輪有效驅動力矩和滾動阻力的時域仿真結果是非穩恒值。經初步分析，認為車輪土壤間的接觸碰撞、車輪輪刺的存在等是產生這一現象的原因。

依據土槽實驗中滑轉率的取值原則，在0~100%選取15個不同值，并通過控制車輪的角速度來控制車輪的滑轉率。根據正交實驗方案，分別進行三個不同車輪荷載W1=100 N，W2=200 N和W3=300 N對應的三組仿真實驗，對所有仿真數據處理后得到圖8(a)所示的三條驅動力矩—滑轉率曲線和圖8(b)所示的三條滾動阻力—滑轉率曲線。為進行仿真結果與實驗結果的對比分析，圖9(a)和(b)繪制了實驗得到的車輪驅動力矩和滾動阻力隨滑轉率的變化曲線。

從實驗結果圖9(a)顯示的驅動力矩—滑轉率曲線可以看出，驅動力矩隨著荷載和(或)滑轉率的增加而增加，相應的仿真結果圖8(a)顯示的驅動力矩隨荷載和滑轉率的變化趨勢與圖9(a)是相同的。同樣，仿真結果圖8(b)顯示的滾動阻力隨荷載和滑轉率的變化趨勢與實驗結果圖9(b)也是一致的。

最后，以200 N荷載下的滾動阻力為例進行了仿真和實驗結果的對比，如圖10所示，圖中用△、符號標志出了實驗結果95%的置信區間。由圖10可以看出，除車輪滑轉率較小(s<0.05)時仿真計算結果具有較小誤差之外，大部分仿真結果具有良好的可信性。

6 結束語

本文基于計算多體系統大規模接觸碰撞動力學理論與地面力學理論相結合的方法，開展了車輪土壤相互作用課題的仿真研究。對仿真結果和實驗結果進行的對比與分析顯示了本仿真工作可以很好地仿真車輪土壤相互作用，具有較高的工程實用性。

車輪土壤相互作用課題的仿真研究在越野行駛和目前廣泛開展的月球車月面行駛性能研究工作中具有重要的意義。相對于其他仿真方法，如離散元法、有限元法等，本文使用的大規模接觸碰撞動力學方法具有計算速度快、工程實用性高的特點。而且本文研制的車輪土壤相互作用模塊基于多體系統仿真平臺，可以方便、快速地建立車輛行走機構等復雜機械系統的多體仿真模型，適合于開展越野車輛、月球車等行駛性能的仿真研究。

參考文獻:

[1]BEKKER M G. Introduction of terrain-vehicle systems[M]. Michigan: The University of Michigan Press， 1969.

[2]WONG J Y. Theory of ground vehicles[M].3rd ed. New York:Wiley， 2001.

[3]莊繼德. 計算汽車地面力學[M]. 北京:機械工業出版社，2002.

[4]ANDRADE G， BEN A F， BIDAUD P， et al. Modeling robot-soil interaction for planetary rover motion control[C]//Proc of IEEE/RSJ International Conference on Intelligent Robots and Systems. 1998: 576-581.

[5]WANJII S， HIROMA T， OTA Y， et al. Prediction of wheel perfor-mance by analysis of normal and tangential stress distributions under the wheel-soil interface[J]. Journal of Terramechanics， 1997， 34(3):165-186.

[6]REINA G， OJEDA L， ANNALISA M， et al. Wheel slippage and sinkage detection for planetary rovers[J]. Trans on Mechatronics， 2006， 11(2):185-195.

[7]NAKASHIMA H， OIDA A. Algorithm and implementation of soil-tire contact analysis code based on dynamic FE-DE method[J]. Journal of Terramechanics， 2004， 41(2-3):127-137.

[8]GODBOLE R， ALCOCK R， HETTIARATCHI D. The prediction of tractive performance on soil surfaces[J]. Journal of Terramecha-nics， 1993， 30(6):443-459.

[9]BLESSING M， DONG S， SMITHBAUER D， et al. Physically based modeling simulation testbed for robotic exploration of a martian surface[C]//Proc of the 1st Space Exploration Conference on Continuing the Voyage of Discovery. 2005.

[10]MAO S G， HAN R P S. Nonlinear complementarity equations for modeling tire-soil interaction:an incremental Bekker approach[J]. Journal of Sound and Vibration， 2006， 312(3):380-398.

[11]GRAHN M. Prediction of sinkage and rolling resistance for off-the-road vehicle considering penetration velocity[J]. Journal of Terramechanics， 1991， 28(4):339-347.

[12]GIBBESCH A， SCHAEFER B， MICHAUD S. 3D simulation and vali-dation of RCL-E and MER rover types mobility[C]//Proc of the 9th ESA Workshop on Advanced Space Technologies for Robotics and Automation. 2006:28-30.

[13]MURILO G. Dynamic simulations of multibody systems[M]. New York: Springer， 2001.

[14]JANOSI Z， HANAMOTO B. Analytical determination of drawbar pull as a function of slip for tracked vehicles in deformable soils[C]//Proc of the 1st International Conference on Terrain-Vehicle Systems. 1961.

[15]MATLAB/Simulink module AESCO soft soil tyre model (ASSTM)user’s guide[R].[S.l.]:AESCO，2003.

[16]孫剛，高峰，李雯. 地面力學及其在行星探測研究中的應用[J]. 力學進展， 2007，37(3):453-464.