車輪-地面耦合動力學行為的CDEM-DEM 分析

路 揚，馮 春，程鵬達，張一鳴

（1.河北工業大學 土木與交通學院，天津 300401；2.中國科學院 力學研究所，北京 100190；3.中國科學院大學 工程科學學院，北京 100049）

1 引言

越野車具有車身為非承載式、車身底盤高、輪胎抓地性好等優點，能夠適用于各種惡劣的道路環境。車輛在軟土路面行駛的過程中，其車輪下陷深度、車輛通行速度、輪胎胎面花紋、地面特征及其相互作用，都會對汽車正常行駛產生一定影響。近年來隨著社會的發展，車輪-地面作用的研究也尤為重要。

針對此類問題，國內外學者進行了大量研究工作。臧孟炎等 人［1-3］基于有限 元（Finite Element Method，FEM）和離散元（Discrete Element Method，DEM）進行了大量仿真分析，通過探究其制動性能和動態特性等，證明了采用數值模擬方法對車輪-地面動態力學行為進行分析的有效性，并對于越野輪胎進行溝槽設計的必要性提出了合理解釋；Zeng 等人［4］基于DEM-FEM 耦合模型通過三軸試驗調整模型參數，證明了該模型對于輪胎-地形相互作用仿真模擬的可行性；Farhadi 等人［5］通過建立輪胎-土壤相互作用的有限元模型并研究了垂直載荷、輪胎充氣壓力、土壤含水量對輪胎功率損失的影響，由此證實了該功率損失估算模型的可行性；鄧露等人［6］基于車輛輪胎相關特征并結合理論推導，提出了一種車輪輪胎的精細化模型，使得其計算結果與實驗數據更加接近準確。

車輪-地面力學動態行為的研究在運輸業、農業、航天等方面均有重要意義，近年來人們對于該方面的研究越發深入。謝斌等人［7］對國內外研究工作進行了綜述，重點討論了輪胎部分的建模方法并結合當今最新研究成果對未來提出了展望；呂鳳天等人［8］通過對車輪-地面進行圖像提取，基于視覺即可對月球車車輪滑轉率進行估計；姜春霞等人［9］基于單輪土槽實驗，對人字形花紋輪胎-地面中的垂直應力分布規律進行了探究；張銳等人［10］基于現場試驗和數值模擬相結合，對車輪曲率半徑在車輪-沙地作用中沉陷性能的影響進行了分析，為沙地車輪設計提供了理論依據。Guo 等人［11］基于多球DE-FE 方法對于某越野氣動輪胎在不規則碎石地形上的牽引性能進行了數值模擬，證明了該方法能夠較好地再現輪胎行駛行為。

綜合以上結果可以看出，數值模擬方法對車輪通行的仿真模擬計算的可行性具有較為完備的理論依據。但由于車輪與軟土路面之間具有較為復雜的動態力學性質，采用單一的有限元方法或離散元方法均無法進行準確刻畫。為此，本文結合前人研究成果，提出了一種CDEM-DEM耦合方法，實現了車輪與軟土路面的動態耦合，并重點探討了車輪花紋情況、路障對車輪行進效率、行進能耗的影響。

2 計算方法介紹

2.1 CDEM 方法

連 續-非 連續單元法［12］（Continuum Discontinuum Element Method，CDEM）是一種有限元與離散元相結合的顯式動力學數值分析方法，其理論基礎為拉格朗日方程，見式（1）：

其中：ui，vi為廣義坐標；L為拉格朗日系統的能量；Qi為非保守力做的功。

連續-非連續單元法中的核心控制方程見式（2）：

其中：M，C，K，Kc，Cc和F分別為單元質量矩陣，單元阻尼矩陣，單元剛度矩陣，接觸面剛度矩陣，接觸面阻尼矩陣和節點外部荷載列陣；分別為單元內所有節點的加速度列陣，速度列陣，位移列陣以及虛擬裂縫上的相對位移列陣，相對速度列陣。

2.2 DEM 方法

DEM［13］是一種顯式數值計算方法，通常用于計算顆粒流在相關條件下的力學性能。DEM的基本原理有兩個方面：接觸模型和運動方程。兩種基本原理分別用來求解離散單元的接觸力和運動狀態物理量（速度、加速度、位移轉角等）。

接觸模型，即力-位移關系。該關系通過離散單元間的連接模型來表征材料本構關系，其理論公式如式（3）、（4）：

運動方程基于牛頓第二定律，其核心控制方程如式（5）、（6）：

其中：mi為剛體質量；rij為j作用于i上的作用點到i形心的距離；vi，ωi分別為i的速度矢量和角速度矢量；Ii為i的慣性矩；g代表重力加速度；t代表時間。

2.3 CDEM-DEM 耦合方法

CDEM 與DEM 分別可對塊體有限元和顆粒離散元的力學行為進行計算，將兩種方法有機耦合，即可對車輪-地面力學動態行為進行準確刻畫。如圖1 為CDEM-DEM 耦合方法示意圖。

圖1 CDEM-DEM 耦合方法Fig.1 CDEM-DEM coupling method

其中，上半部分CDEM 模型包括塊體和界面兩部分，塊體由有限元單元體所組成；界面為塊體間的公共邊界，分別表征材料的彈性、塑性、損傷等連續特征和斷裂、滑移、碰撞等非連續特征。數值模擬時，車輪采用CDEM 塊體進行描述，軟土路面采用DEM 顆粒進行描述。塊體與顆粒、顆粒與剛性面之間接觸耦合剛度采用全局的值，法向、切向剛度均為1×109Pa·m-1。節點力和節點運動的核心控制方程見式（7）、（8）。基于以下方程即可實現顯式求解過程。

其中：F為節點合力，FE為節點外力，Fe為有限元單元變形貢獻的節點力，Fc為接觸界面貢獻的節點力，Fd為節點阻尼力，a為節點加速度，v為節點速度，Δu為節點位移增量，u為節點位移全量，m為節點質量，Δt為計算時步。

CDEM-DEM 耦合方法主要包含接觸檢測與接觸力計算兩部分，其中接觸力計算與顆粒離散元間接觸力［14］一致，均為基于增量法的顯式求解方式，其計算式見式（9）、（10）：

其中：Fn，Fs分別為顆粒間的法向和切向接觸力；Δt為計算時步；Δdun，Δdus分別為兩個接觸顆粒間的法向和切向位移增量差。

接觸檢測首先基于子空間法，通過顆粒、塊體等之間的映射關系縮小檢測范圍，之后，進行精確檢測。在本文的車輪-地面問題中，該接觸檢測即是有限元與邊界面之間的接觸，一般采用點-面接觸模型［15］實現。其需要滿足如下兩個基本條件：

其中：d為顆粒體心到邊界面的距離，R為顆粒球體半徑，C為邊界面，O為顆粒體心在邊界面上的投影位置點。

2.4 動態邊界條件施加

由于CDEM 中的單元只有平動自由度，沒有轉動自由度，無法直接對CDEM 單元施加扭矩；為此，本文提出了一種在局部坐標系下施加動態切應力，從而間接施加動態扭矩的方法。在車輪輪轂面施加線性變化的動態邊界條件從而牽引車輪轉動，開啟坐標系施加開關，在兩個切向上施加動態邊界條件。其核心控制方程如式（13）：

其中：Fi為第i方向的載荷值（面力）；fT0，fT1分別為線段起始時間和結束時間；fV0，fV1分別為線段起始值和線段結束值；λi為各個方向的載荷系數；t為當前時間。

圖2 為動態邊界條件施加方式示意圖，在該模型中，通過在車輪輪轂中心圓柱面設定坐標系以確定法向與切向，在圓柱面上施加切向面力從而產生轉矩使車輪轉動。

圖2 動態邊界條件施加方式示意圖Fig.2 Schematic diagram of application mode of dynamic boundary conditions

通過在圓柱面每一個節點施加面力，對所有面力進行求和形成轉矩，其主要控制方程如式（14）、（15）：

其中：Fi為作用面力，Ai為面力作用的面積，τ為切應力，M為轉矩，R為圓柱面半徑。

轉速力矩控制方程如式（16）、（17），通過推導得出轉速切應力關系如式（18）：

其中：β為角加速度，J為轉動慣量，m為質量，ω為轉速，t為時間，ρ為密度。

創建圓柱體車輪進行驗證，車輪截面半徑R為0.5 m，高h為 0.2 m，密度為1 000 kg·m-3，在圓柱體外環施加局部坐標系下動態切應力1 000 Pa，根據上述方程求得轉速預測結果。圖3 為預測與模擬結果對比圖。

圖3 預測與模擬結果對比Fig.3 Comparison of prediction and simulation results

數值計算結果與理論解基本一致，表明了本文提出的利用局部坐標系下切向應力施加扭矩的計算方法的正確性。

3 仿真模擬分析

3.1 模型建立及材料參數

車輪模型尺寸采用了8.25R16LT 132/128F 18PR 全鋼載重子午線輪胎［16］。該輪胎外直徑855 mm，斷面寬230 mm，為簡化計算，同時使車輪嵌入地面時車轍更加明顯，花紋節距數采用20，其余參數不變。花紋車輪模型見圖4，同時建立如圖5 所示光面車輪進行對比分析。

圖4 花紋車輪模型Fig.4 Patterned wheel model

圖5 光面車輪模型Fig.5 Glossy wheel model

其中，輪胎材料選取隔振橡膠［17］，采用linear線彈性本構。對輪胎材料進行均一化假設，將其彈性模量調整至充氣狀態標準；軟土路面采用石灰 土［18］，采 用brittleMC 脆 性斷裂 本構。石灰土模型長3 m，寬1 m，高0.14 m。材料力學參數見表1。

表1 材料力學參數Tab.1 Material mechanical parameters

假定汽車重量為3.30×103kg，在輪胎中心建立一組同心圓柱輪轂模型以增加車輪重量使之達到四分之一車重，故需進行等效施加較大轉矩使其車速符合實際情況。經實驗，當局部坐標系下動態切應力大小在全程恒為10 MPa 時，車輪轉速將穩定于11.00 m·s-1，適用于汽車在軟土路面行駛時的車速。

在此基礎上，在模型周圍創建5 個剛性面進行約束。為便于計量，記花紋輪胎工況為T1，光面輪胎工況記為T2，二者初始埋深約為0.03 m。其最終計算模型分別如圖6 和圖7 所示。

圖6 T1 計算模型Fig.6 T1 calculation model

圖7 T2 計算模型Fig.7 T2 calculation model

3.2 模型驗證

如圖8 為本次模擬結果與實驗結果對比圖（以T2 為例），在計算過程中，與車輪接觸的軟土路面分為如圖8（b）中兩個區域，即前區順時針區和后區逆時針區。如圖，T2 在t=0.40 s 時產生顆粒飛濺現象，其前方軟土運動方向為順時針；后方軟土產生逆時針方向運動趨勢導致車輪后方軟土產生堆積現象。該現象與文獻［2］中所述實驗結果相同。

圖8 模擬結果與實驗結果對比Fig.8 Comparison between simulation result and experimental result

由此可見該計算方法可準確模擬車輛在軟土路面行駛的真實情況。

4 計算結果

4.1 T1、T2 工況對比

本次計算過程中，車輪在石灰土軟土路面上行駛，待車輪轉動至終點計算完畢。如圖9～10為計算的兩組工況車轍圖。從圖中可以看到，軟土路面由于車輪駛過產生顆粒的擠壓流動導致車輪兩側和后側產生顆粒隆起現象，在地面上留下了清晰可見的車轍。在相同時間內，T1 在軟土路面運行完畢到達終點共耗時1.90 s，其運行長度為2.15 m；而T2 在此期間僅運行了0.72 m。由此可見在軟土路面上，花紋輪胎的通行能力較光面輪胎更為出色，相同時間內其通行路程更遠。

圖9 T1 車轍Fig.9 T1 rutting

圖10 T2 車轍Fig.10 T2 rutting

圖11 為車輪速度時程圖，從圖中可以看出，在施加動態切應力后，T1、T2 中車輪開始轉動，其轉動速度與平動速度逐漸升高，而后趨于穩定。

圖11 車輪速度時程圖Fig.11 Wheel speed-time diagram

其中，T1 轉動速度穩定值約為21.05 rad·s-1，T2 轉動速度穩定值約為23.39 rad·s-1，其對應轉動線速度分別為9.00 m·s-1，10.00 m·s-1；T1 平動速度穩定值約為1.15 m·s-1，T2 平動速度穩定值約為0.28 m·s-1。

光面車輪由于其表面光滑，與地面滑動摩擦力較小，故在轉動過程中，其轉動速度略高于花紋車輪，車輪與地面相對運動時發生打滑現象，故其通行能力遠低于花紋車輪。在T1 中，花紋車輪平動速度與其線速度之比約為12.78%；T2中光面車輪的比值約為2.80%。根據該比值可看出，在軟土路面行駛過程中，相較于光面車輪，花紋車輪通行能力更好，由此說明了在汽車制造設計中輪胎胎面合理設計的必要性。

如圖12 為系統動能時程圖，從圖中可以看到T1 與T2 中整個系統動能在0～0.50 s 內均持續上升，但T2 動能增長速度明顯高于T1，這是因為T2 中的光面輪胎由于其與地面滑動摩擦力小，導致其轉動速度高于T1 中的花紋車輪。故在此期間T2 動能增長速度高于T1。在0.50 s之后，T1、T2 系統動能均趨于穩定，其穩定值分別約為6 000 J、10 000 J。之所以出現這種差異是因為在此過程中，光面車輪相較于花紋車輪，其下陷深度更大，故在此過程中，車輪轉動所引起的顆粒流動現象更為劇烈，運動顆粒范圍較T1 更大，因此T2 系統動能大于T1。由此可見，車輪下陷程度也是影響車輛在軟土路面通行能力的重要因素之一。

圖12 動能時程圖Fig.12 Kinetic energy time history

圖13 和 圖14 分別為T1、T2 在Y方向的位移云圖。從圖中可以看出，兩組工況中由于車輪前進，在車輪后方產生了顆粒向后流動的現象：在t=1.00 s 時，T1 車輪后產生顆粒后移現象，其位移大小約為0.70 m；T2 約為0.40 m。

圖13 T1Y方向位移云圖Fig.13 Displacement cloud diagram of T1 onYdirection

圖14 T2Y方向位移云圖Fig.14 Displacement cloud diagram of T2 onYdirection

在t=1.90 s 計算完畢時，T1 位移大小約為0.85 m；T2 約為0.55 m。不難看出，T1 中在花紋車輪作用下顆粒后移數值較大，在圖中也可看出在T1 中發生了較為明顯的顆粒飛濺現象（軟土路面行駛出現塵土飛揚）。由此也可說明在軟土路面中，花紋車輪相較于光面車輪轉動會產生較為明顯的軟土后移現象從而產生軟土堆積，驅動車輪向前運動，由此提高了車輪通行能力。

在圖9 和圖10 中可以看到，隨著車輪在軟土路面的運行，車輪下陷深度出現明顯差異，花紋車輪下陷深度低于光面車輪，如圖15 為車輪下陷深度時程圖。

圖15 車輪下陷深度時程圖Fig.15 Time history of wheel sinking depth

從圖中可以看出，T1 中花紋車輪隨著時間的推進，在0～0.50 s 中，其下陷深度先升高后降低，之后在0.50 s 時車輪下陷深度約為0.01 m，之后便穩定于該值；T2 中光面車輪在0～0.75 s內下陷深度不斷增大，在達到0.12 m 后穩定于該值。由此可見，在質量相同的情況下，花紋車輪的下陷深度遠低于光面車輪。故合理地進行車輪輪胎胎面設計，可以減小車輪侵入軟土路面的深度，提高通行能力。

4.2 含路障路面探究

為探究車輪在含路障路面的通行能力，設置三組路障，路障高度分別為0.02 m，0.05 m，0.07 m，記三種工況分別為U1、U2、U3，圖16 為三種路障示意圖。

圖16 路障示意圖Fig.16 Schematic diagram of roadblock

花紋車輪在三種工況中行駛1 s，最終在U1、U2、U3 的行駛距離分別為1.04 m，0.96 m，0.95 m。U1 工況通行距離高于后兩種工況，U2、U3通行距離相近。

圖17 為在路障工況下的系統動能時程圖，從圖中可以看出隨著時間的推移，三種工況動能逐漸上升。

圖17 路障工況下動能時程圖Fig.17 Kinetic energy time history under roadblock condition

與T1、T2 工況中動能隨著時間的推移趨于穩定不同，U1、U2、U3 在0.5 s 后動能值出現明顯波動。三種工況動能均出現先升高后下降的趨勢：U1、U2、U3 均大約在0.8 s 時達到峰值，此時U1、U2、U3 的動能分別 為7 815.11 J、9341.75 J、9 822.38 J。說明在此過程中，車輪遇到障礙，必須增大制動力才能通過，由此動能上升，在通過障礙物后，動能逐漸下降。軟土障礙越大，其所需能量越大，持續時間越長。

由此可見車輛在含路障路面行駛的過程中，面對障礙物需要增強制動力以通過障礙物，也要在通過障礙物后適當減速以保持車輛平穩通過。

5 力學機理探究

綜合以上計算結果可知，花紋輪胎由于其胎面表面凹凸不平，增大了與軟土路面的接觸面積，從而增大了滑動摩擦力，故花紋車輪轉速低于光面車輪。但由于胎面花紋的存在，使花紋車輪產生切向作用力，在一定程度上增大了花紋車輪的驅動力，故花紋車輪不會發生如光面車輪“打滑下陷”的現象，而是保持勻速前進。光面車輪由于驅動力不足無法順利通行，其平動速度低于花紋車輪，但其轉動速度較高，與顆粒發生較為劇烈的相對運動，從而導致軟土路面松軟，故其下陷深度遠高于花紋車輪。相同時間內，花紋車輪通行路程更遠，通行能力更強；車輪在含路障路面行駛時，往往會有坡度高低不同的路障，在車輪通過路面障礙時，車輛需要增強驅動力以獲取足夠能量跨越障礙物，尤其是起伏程度較大的路障，其所需能量往往更大。在車輪通過障礙物后其能耗又會有所下降，這就需要車輛在行駛過程中進行合理制動從而平穩通過。

6 結論

為了準確刻畫車輪在軟土路面上轉動、摩擦、前行的動力學行為，提出了CDEM-DEM 耦合的數值計算方法，實現了通過施加局部坐標系下動態切向面力間接實現動態轉矩施加的方法，探討了車輪花紋、路障對車輪前行速率及能耗的影響規律。研究結果表明：

（1）花紋輪胎、光面輪胎均在軟土路面留下了清晰可見的車轍；相同質量下，花紋車輪下陷深度遠低于光面車輪下陷深度。

（2）花紋車輪與光面車輪在轉動前行時均會產生較為明顯的軟土后移現象，其軟土后移值分別為0.85 m、0.55 m。

（3）在相同扭矩情況下，光面車輪轉速更快，并出現了打滑前行的現象；相反，花紋車輪的轉速較慢，但在軟土路面上卻以較高的速度前行；花紋車輪平動速度與其線速度之比約為12.78%；而光面車輪比值僅為2.80%。

（4）輪胎在含路障路面通行時，其地面起伏程度越大，所消耗的能量越大。這就要求汽車在軟土路面行駛時，應對于路況的不同做出合理的制動以平穩通過。