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

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

1 引言

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

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

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

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

2 計算方法介紹

2.1 CDEM 方法

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

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

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

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

2.2 DEM 方法

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

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

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

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

2.3 CDEM-DEM 耦合方法

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

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

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

其中：F為節(jié)點合力，F(xiàn)E為節(jié)點外力，F(xiàn)e為有限元單元變形貢獻的節(jié)點力，F(xiàn)c為接觸界面貢獻的節(jié)點力，F(xiàn)d為節(jié)點阻尼力，a為節(jié)點加速度，v為節(jié)點速度，Δu為節(jié)點位移增量，u為節(jié)點位移全量，m為節(jié)點質(zhì)量，Δt為計算時步。

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

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

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

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

2.4 動態(tài)邊界條件施加

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

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

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

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

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

其中：Fi為作用面力，Ai為面力作用的面積，τ為切應(yīng)力，M為轉(zhuǎn)矩，R為圓柱面半徑。

轉(zhuǎn)速力矩控制方程如式（16）、（17），通過推導(dǎo)得出轉(zhuǎn)速切應(yīng)力關(guān)系如式（18）：

其中：β為角加速度，J為轉(zhuǎn)動慣量，m為質(zhì)量，ω為轉(zhuǎn)速，t為時間，ρ為密度。

創(chuàng)建圓柱體車輪進行驗證，車輪截面半徑R為0.5 m，高h為 0.2 m，密度為1 000 kg·m-3，在圓柱體外環(huán)施加局部坐標系下動態(tài)切應(yīng)力1 000 Pa，根據(jù)上述方程求得轉(zhuǎn)速預(yù)測結(jié)果。圖3 為預(yù)測與模擬結(jié)果對比圖。

圖3 預(yù)測與模擬結(jié)果對比Fig.3 Comparison of prediction and simulation results

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

3 仿真模擬分析

3.1 模型建立及材料參數(shù)

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

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

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

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

表1 材料力學(xué)參數(shù)Tab.1 Material mechanical parameters

假定汽車重量為3.30×103kg，在輪胎中心建立一組同心圓柱輪轂?zāi)Ｐ鸵栽黾榆囕喼亓渴怪_到四分之一車重，故需進行等效施加較大轉(zhuǎn)矩使其車速符合實際情況。經(jīng)實驗，當(dāng)局部坐標系下動態(tài)切應(yīng)力大小在全程恒為10 MPa 時，車輪轉(zhuǎn)速將穩(wěn)定于11.00 m·s-1，適用于汽車在軟土路面行駛時的車速。

在此基礎(chǔ)上，在模型周圍創(chuàng)建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 為本次模擬結(jié)果與實驗結(jié)果對比圖（以T2 為例），在計算過程中，與車輪接觸的軟土路面分為如圖8（b）中兩個區(qū)域，即前區(qū)順時針區(qū)和后區(qū)逆時針區(qū)。如圖，T2 在t=0.40 s 時產(chǎn)生顆粒飛濺現(xiàn)象，其前方軟土運動方向為順時針；后方軟土產(chǎn)生逆時針方向運動趨勢導(dǎo)致車輪后方軟土產(chǎn)生堆積現(xiàn)象。該現(xiàn)象與文獻［2］中所述實驗結(jié)果相同。

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

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

4 計算結(jié)果

4.1 T1、T2 工況對比

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

圖9 T1 車轍Fig.9 T1 rutting

圖10 T2 車轍Fig.10 T2 rutting

圖11 為車輪速度時程圖，從圖中可以看出，在施加動態(tài)切應(yīng)力后，T1、T2 中車輪開始轉(zhuǎn)動，其轉(zhuǎn)動速度與平動速度逐漸升高，而后趨于穩(wěn)定。

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

其中，T1 轉(zhuǎn)動速度穩(wěn)定值約為21.05 rad·s-1，T2 轉(zhuǎn)動速度穩(wěn)定值約為23.39 rad·s-1，其對應(yīng)轉(zhuǎn)動線速度分別為9.00 m·s-1，10.00 m·s-1；T1 平動速度穩(wěn)定值約為1.15 m·s-1，T2 平動速度穩(wěn)定值約為0.28 m·s-1。

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

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

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

圖13 和 圖14 分別為T1、T2 在Y方向的位移云圖。從圖中可以看出，兩組工況中由于車輪前進，在車輪后方產(chǎn)生了顆粒向后流動的現(xiàn)象：在t=1.00 s 時，T1 車輪后產(chǎn)生顆粒后移現(xiàn)象，其位移大小約為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 中在花紋車輪作用下顆粒后移數(shù)值較大，在圖中也可看出在T1 中發(fā)生了較為明顯的顆粒飛濺現(xiàn)象（軟土路面行駛出現(xiàn)塵土飛揚）。由此也可說明在軟土路面中，花紋車輪相較于光面車輪轉(zhuǎn)動會產(chǎn)生較為明顯的軟土后移現(xiàn)象從而產(chǎn)生軟土堆積，驅(qū)動車輪向前運動，由此提高了車輪通行能力。

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

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

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

4.2 含路障路面探究

為探究車輪在含路障路面的通行能力，設(shè)置三組路障，路障高度分別為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 為在路障工況下的系統(tǒng)動能時程圖，從圖中可以看出隨著時間的推移，三種工況動能逐漸上升。

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

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

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

5 力學(xué)機理探究

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

6 結(jié)論

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

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

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

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

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