程鵬達(dá)，馮 春，2*，常寧東，周 俊，朱心廣

（1.中國科學(xué)院 力學(xué)研究所，北京 100190；2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院，北京 100049）

1 引言

車輛通行性研究在軍事、建筑業(yè)和農(nóng)業(yè)等領(lǐng)域起著重要作用，尤其在軍事方面，高機(jī)動性和可靠性已成為現(xiàn)代車輛設(shè)計和研究的重要目標(biāo)。當(dāng)場地環(huán)境復(fù)雜時，車輛通行性不僅與車輛本身有關(guān)，還與地面地形條件和地面土力學(xué)特性相關(guān)。因此，車輛-地形力學(xué)耦合動態(tài)特性的研究具有重要的應(yīng)用價值和科學(xué)意義。

車輛車身（如輪胎、履帶、板、腿、底座、葉片和尖齒）與土體、土壤嵌入障礙物（如植被根/莖、巖石）接觸過程對車輛通行性有重要影響。近年來，美國陸軍、北約成員國和伙伴國的軍隊持續(xù)提出車輛全地形通行性預(yù)測的需求［1-2］。當(dāng)前，北約采用的車輛通行性評估主要側(cè)重于車輛機(jī)動性，即車輛速度（NATO Reference Mobility Model，NRMM）模型［3］。NRMM 模型開發(fā)于20 世紀(jì)70 年代，并在80 年代和90 年代持續(xù)更新，該模型基于經(jīng)驗關(guān)系來預(yù)測給定地形上的車輛機(jī)動性，同時考慮地形力學(xué)等變量，如土壤圓錐指數(shù)（CI）［4］、USCS（ASTM 2000）土壤類型、土壤濕度、地表覆蓋條件（正常、水或雪）和深度、土地利用（農(nóng)業(yè)、灌木叢、森林、濕地、城市等）、坡度（上坡、下坡和側(cè)坡）、表面粗糙度、土堆/溝渠障礙物尺寸和間距、樹木/植被莖的大小和間距、可見度以及道路參數(shù)（道路類型、車道、橋梁、隧道等）?；诮?jīng)驗關(guān)系，NRMM 對20 世紀(jì)60 年代至80 年代的車輛適用性較好，但很難推廣到當(dāng)前多類型和尺寸的車輛，同時調(diào)整經(jīng)驗關(guān)系需要全尺寸車輛試驗，非常昂貴且耗時。21 世紀(jì)初，北約提出需要更新、改進(jìn)和重新驗證軟土-車輛耦合動力學(xué)模型。2016 年至2018 年，北約RTG-248 研究任務(wù)組提出分階段更新方案，短期方案基于Bekker-Wong 理論［5-6］，該理論較Mohr-Coulomb 模型能更好地描述土壤剪切應(yīng)力應(yīng)變關(guān)系，也稱之為簡單地形力學(xué)模型。長期方案是基于離散元法（Discrete Element Method，DEM）或有限元（Finite Element Method，F(xiàn)EM）模型，也稱之為復(fù)雜地形力學(xué)模型。基于Bekker-Wong 理論，國內(nèi)外眾多學(xué)者開展了大量研究［7-10］，針對不同類型車輛-土壤相互作用簡化關(guān)系，推導(dǎo)了二維或三維條件下土壤特性、剪切速率、土壤推力、牽引力、阻力等參數(shù)間關(guān)系?；趶?fù)雜地形力學(xué)模型，能夠更為真實地預(yù)測車輛車身（輪胎、履帶、支腿、鏟斗、葉片、尖齒等）上的三維反作用力，以及土壤三維變形和流動［1］，如車轍深度/寬度/形狀、車轍側(cè)壁高度和車輛前方的推土。與車輛耦合的復(fù)雜地形力學(xué)通?；谶B續(xù)介質(zhì)模型、連續(xù)-非連續(xù)模型、非連續(xù)模型，常用方法包括有限元法（FEM）［11-12］、連續(xù)-非連續(xù)單元法（Continuum-discontinuum Element Method，CDEM）、離散元法（DEM）、光滑粒子流體動力學(xué)（Smooth Particle Hydrodynamics，SPH）［13］和物質(zhì)點(diǎn)法（Material Point Method，MPM）等。車輛-土壤耦合作用的物理現(xiàn)象十分復(fù)雜，通常借助數(shù)值方法，以較高的計算成本為代價得到更為詳細(xì)的結(jié)果，且不依賴經(jīng)驗參數(shù)。值得注意的是粒子模型將土壤模擬為粒子集合，可以捕捉到土壤變形流動的復(fù)雜過程［11］，例如推土或土壤混合過程，而連續(xù)介質(zhì)模型對該過程的描述存在一定困難。此外，針對土壤應(yīng)力應(yīng)變本構(gòu)關(guān)系，Bekker 做出了開創(chuàng)性研究［5］，Wong 等在Bekker 本構(gòu)關(guān)系基礎(chǔ)上提出了Bekker-Wong 本構(gòu)關(guān)系［6］。近年來，許多研究人員專注于更精確的半經(jīng)驗本構(gòu)關(guān)系，多數(shù)情況下這些本構(gòu)關(guān)系仍然受到Bekker 模型的啟發(fā)［14］。

目前關(guān)于車輛-地形力學(xué)多介質(zhì)耦合動力學(xué)仍是車輛裝備通行性研究的核心問題，具有重要的應(yīng)用價值和科學(xué)意義。因此，本文基于自主知識產(chǎn)權(quán)的連續(xù)-非連續(xù)單元法（CDEM）［15-16］和物質(zhì)點(diǎn)法（MPM），數(shù)值模擬某型裝備引起軟土地面的變形及侵入軟土的過程，定量研究車輪載重、軟土彈性模量、軟土強(qiáng)度參數(shù)與軟土表面壓應(yīng)力、侵入深度等參數(shù)之間的關(guān)系，獲得車輪侵入軟土深度的變化規(guī)律，為復(fù)雜地形土壤環(huán)境下軍事行動分析、任務(wù)規(guī)劃和戰(zhàn)爭模擬預(yù)測的車輛通行性提供理論依據(jù)。

2 數(shù)值方法

2.1 連續(xù)-非連續(xù)單元法（CDEM）

基于拉格朗日方程的連續(xù)-非連續(xù)單元法（CDEM）是一種有限元與離散元耦合的顯式數(shù)值分析方法，主要用于材料漸進(jìn)破壞過程的模擬。CDEM 的計算域通常包含連續(xù)塊體和分散塊體，分別作為有限元（FEM）和離散元（DEM）方法的計算域［15-16］，如圖1 所示。

圖1 CDEM 計算域Fig.1 Computational domain in the CDEM

CDEM 基于增量方式的顯式歐拉前差法求解動力問題，主要包含節(jié)點(diǎn)合力計算及節(jié)點(diǎn)運(yùn)動計算兩個部分。節(jié)點(diǎn)合力如下：

其中：F為節(jié)點(diǎn)合力，F(xiàn)E為節(jié)點(diǎn)外力，F(xiàn)e為節(jié)點(diǎn)變形力（由單元應(yīng)力貢獻(xiàn)），F(xiàn)c為接觸界面貢獻(xiàn)的節(jié)點(diǎn)力，F(xiàn)d為節(jié)點(diǎn)阻尼力。

節(jié)點(diǎn)運(yùn)動計算如下：

其中：a為節(jié)點(diǎn)加速度，v為節(jié)點(diǎn)速度，u為節(jié)點(diǎn)位移全量，m為節(jié)點(diǎn)質(zhì)量，Δt為計算時步。

CDEM 中，連續(xù)求解域的有限元單元應(yīng)力及節(jié)點(diǎn)變形力采用增量法計算：

其中：Bi，Δεi，Δσi，wi和Ji分別為高斯點(diǎn)i的應(yīng)變矩陣，增量應(yīng)變向量，增量應(yīng)力向量，積分系數(shù)及雅克比行列式；n和n-1 為高斯點(diǎn)i當(dāng)前時刻及上一時刻；D，Δue，F(xiàn)e分別表示單元的彈性矩陣，節(jié)點(diǎn)增量位移及節(jié)點(diǎn)力；N表示高斯點(diǎn)個數(shù)。

CDEM 中，將非連續(xù)求解域劃分為不同于原始連續(xù)塊體的分散塊體，這些分散塊體之間的相互作用可以轉(zhuǎn)化為虛擬彈簧接觸力。本文采用半彈簧-半棱聯(lián)合接觸模型快速標(biāo)記接觸對，并計算接觸力［15-18］。

局部坐標(biāo)系下，塊體間的接觸力可以表達(dá)為：

其中：F，K和Δu分別為虛擬彈簧力的增量，剛度和相對位移；下標(biāo)nt 和s 分別代表法向和切向。

基于不同介質(zhì)的物理力學(xué)特性，當(dāng)虛擬彈簧達(dá)到給定的失效準(zhǔn)則時，塊體間連續(xù)界面將轉(zhuǎn)化為不連續(xù)的斷裂面，用于表征材料的斷裂、滑移、碰撞等非連續(xù)特征。本文選擇最大拉應(yīng)力準(zhǔn)則為失效準(zhǔn)則，剪切強(qiáng)度選擇Bekker 模型描述，對于準(zhǔn)靜態(tài)問題，該模型可以簡化為Mohr-Coulomb 模型。

2.2 物質(zhì)點(diǎn)法（MPM）

物質(zhì)點(diǎn)法（MPM）是由Sulsky 等人提出的一種無網(wǎng)格的方法［19］。MPM 方法將連續(xù)體離散成一組質(zhì)點(diǎn)，質(zhì)點(diǎn)攜帶了所有物質(zhì)信息，其運(yùn)動代表物質(zhì)的變形，避免了界面描述等問題。質(zhì)點(diǎn)和背景網(wǎng)格固連后，物質(zhì)點(diǎn)信息映射到背景網(wǎng)格點(diǎn)，動量方程在背景網(wǎng)格上求解，得到網(wǎng)格點(diǎn)加速度、速度和位移。將網(wǎng)格點(diǎn)運(yùn)動量映射回物質(zhì)點(diǎn)，得到物質(zhì)點(diǎn)的運(yùn)動信息。丟棄變形后的背景網(wǎng)格，重新構(gòu)造背景網(wǎng)格，進(jìn)而避免了網(wǎng)格畸變。MPM 方法結(jié)合Lagrange 方法和Euler 方法的優(yōu)點(diǎn)，適用于沖擊、碰撞、破壞、加工成型等大變形問題，但是反復(fù)映射帶來較大的計算量，并容易產(chǎn)生誤差。

MPM 方法動量方程的弱形式為［20］：

其中：δui、ρ、¨、σij、fi和ti分別為虛位移、密度、加速度、Cauchy 應(yīng)力、體力和邊界牽引力，V表示材料域，Γ表示材料邊界。

離散顆粒的質(zhì)量密度可以表示為：

其中：Mp為質(zhì)量，np為顆?？倲?shù)，δ為狄拉克函數(shù)（Dirac Delta Function），xpi為顆粒p的坐標(biāo)。

將質(zhì)量密度代入動量方程的弱形式，得到：

其中，h-1是滿足體積積分所需的邊界層假設(shè)厚度［21］。

MPM 方法的求解過程中，質(zhì)點(diǎn)和背景網(wǎng)格固連。每個時間步計算結(jié)束后，更新粒子的速度和位置，并為下一個時間步定義新的規(guī)則網(wǎng)格。因此，與有限元法相比，避免了網(wǎng)格畸變。通常情況下，所有時間步中都可以使用相同的固定規(guī)則網(wǎng)格。因此，可以通過背景網(wǎng)格中節(jié)點(diǎn)的有限元形函數(shù)，實現(xiàn)粒子和節(jié)點(diǎn)之間的映射：

其中：uip、uiI、uip，j分別為顆粒位移、網(wǎng)格位移和顆粒位移的導(dǎo)數(shù)。g為有限元的網(wǎng)格節(jié)點(diǎn)總數(shù)，對于一個六面體單元，g=8。

2.3 CDEM-MPM 耦合計算方法

本文選擇CDEM-MPM 耦合方法計算車輪和軟土地面相互作用的過程，以精確模擬軟土地面變形、損傷和破壞過程中涉及的多種物理現(xiàn)象。通過引入顆粒-表面/邊緣接觸算法，實現(xiàn)MPM 計算域和CDEM 計算域相互作用，該接觸模型耦合分析如圖2 所示［22］。

圖2 CDEM-MPM 耦合分析示意圖Fig.2 Schematics of CDEM-MPM coupled analysis

如圖2（a）所示，通過滿足以下要求建立接觸：

其中，d和R分別表示距離和顆粒半徑。

當(dāng)CDEM 單元與MPM 物質(zhì)點(diǎn)接觸時，如圖2（b）所示，法向接觸力和剪切接觸力計算分別如下：

其中：k為接觸剛度，n為單元邊界方向，vm為顆粒的速度矢量，vn為單元邊界上顆粒m投影點(diǎn)的速度：

其中：vB和vC為節(jié)點(diǎn)速度。

因此，可以通過顆粒-表面/邊緣接觸法實現(xiàn)CDEM-MPM 的耦合。在每個時間步中，首先檢測MPM 計算域和CDEM 計算域是否接觸，若接觸，則計算它們間的接觸力，實現(xiàn)相互作用；若未接觸，則它們會獨(dú)立更新，以獲取節(jié)點(diǎn)變量的信息。為了在數(shù)值計算中獲得穩(wěn)定的解，背景單元應(yīng)至少包含一個顆粒，即顆粒半徑不超過背景網(wǎng)格單元長度的一半。

CDEM-MPM 耦合方法在解決大變形問題上具有突出優(yōu)勢，同時也大幅降低了計算成本。近年來，相關(guān)研究也驗證了通過顆粒-表面/邊緣接觸法算法實現(xiàn)塊體和粒子耦合（CDEMMPM）算法的準(zhǔn)確性［23-24］。

3 車輪-軟土耦合動力學(xué)數(shù)值模擬

本文采用連續(xù)-非連續(xù)（CDEM）方法對車輪進(jìn)行建模，采用物質(zhì)點(diǎn)法（MPM）對軟土進(jìn)行建模。最后，通過數(shù)值模擬分析車輪-軟土地面的相互作用，并討論車輪侵入深度的變化規(guī)律。在國際地面車輛系統(tǒng)學(xué)會（ISTVS 1977）標(biāo)準(zhǔn)中，車輪侵入深度被定義為從車輪的最低點(diǎn)到未擾動的土壤表面的距離。

3.1 幾何模型

選擇某型工程車常用輪胎，輪胎寬度為52.07 cm，輪胎外徑為180 cm，如圖3 所示。軟土地面為規(guī)則長方體區(qū)域，長度為6 m，寬度為3 m，高度為2 m，如圖4 所示，同時選取軟土表面為零平面。經(jīng)過顆粒尺寸敏感性分析，軟土地面物質(zhì)點(diǎn)顆粒直徑為輪胎寬度的1/32。

圖3 車輪幾何模型Fig.3 Geometry model of wheel

圖4 車輪-軟土侵入計算幾何模型Fig.4 Geometry model of wheel and soft-soil interaction calculation

3.2 物理力學(xué)參數(shù)、初始和邊界條件

選取某黏性土作為典型軟土地面，該黏性土物理力學(xué)參數(shù)如表1 所示［25］。分析不同車輪載重（1 000 kg、3 000 kg、5 000 kg、7 000 kg、9 000 kg和11 000 kg）、彈性模量（4 MPa、8 MPa、12 MPa、16 MPa、20 MPa 和24 MPa）、黏 聚力（5 kPa、10 kPa、15 kPa、20 kPa 和25 kPa）和摩擦角（10°、15°、20°、25°和30°）時，車輪侵入軟土地面的深度。軟土地面力學(xué)參數(shù)變化涵蓋了《工程勘察通用規(guī)范》（GB 55017-2021）中黏性土、粉土和部分砂土的力學(xué)參數(shù)取值范圍［26］。計算過程中，每次變化一個參數(shù)，其余參數(shù)按表1 中數(shù)值取值。當(dāng)車輪氣壓達(dá)標(biāo)時，車輪變形量比軟土變形量低二到三個量級，車輪可采用線彈性本構(gòu)模型。當(dāng)車輪變形遠(yuǎn)小于軟土變形，可簡化車輪中橡膠、鋼絲、簾布以及帶束等材料各自的力學(xué)特性，采用加權(quán)平均后整體車輪的力學(xué)特性，彈性模量大于470 MPa，泊松比約為0.48。

表1 物理力學(xué)參數(shù)Tab.1 Physical and mechanical parameters

軟土地面的初始地應(yīng)力與土壤密度、深度相關(guān)，且初始位移為0。軟土地面四周和底邊界采用法向位移約束，位移為0。

4 計算結(jié)果

4.1 車輪載重的影響

改變車輪載重（1 000 kg、3 000 kg、5 000 kg、7 000 kg、9 000 kg 和11 000 kg），分析車輪-軟土相互作用規(guī)律。以3 000 kg 和11 000 kg 載重為例，車輪侵入軟土的位移云圖分別如圖5 和圖6所示。當(dāng)載重為3 000 kg 時，車輪侵入軟土深度大約為4.61 cm。當(dāng)載重為11 000 kg 時，車輪侵入軟土深度大約為31.27 cm。隨載重增加，車輪侵入深度逐漸增大，車輪徑向（進(jìn)行方向）和軸向軟土受擠壓向上隆起。軟土隆起范圍和高度與載重正相關(guān)，車輪侵入軟土過程中，三維效應(yīng)非常明顯。

圖6 車輪侵入軟土位移云圖（11 000 kg）Fig.6 Displacement of soft-soil under wheel load（11 000 kg）

不同載重條件下，軟土表面壓應(yīng)力在車輪徑向和軸向的差異較大。選取車輪徑向?qū)ΨQ面在軟土表面的投影線及其延伸線為特征線，車輪中心的投影點(diǎn)為0 m 位置。載重變化后，特征線上的軟土表面壓應(yīng)力曲線如圖7 所示。不同載重時，車輪對軟土的壓應(yīng)力均以車輪中心投影點(diǎn)為對稱點(diǎn)接近對稱分布，對稱點(diǎn)（車輪中心的投影點(diǎn)）壓應(yīng)力最大，壓應(yīng)力輪廓隨載重增加而擴(kuò)大，說明車輪接觸軟土的面積逐漸增大。隨載重增加，壓應(yīng)力最大值逐漸增大，分別為15.39 kPa、31.47 kPa、46.09 kPa、59.19 kPa、70.30 kPa 和78.65 kPa。遠(yuǎn)離對稱點(diǎn)時，壓應(yīng)力先快速降低后小幅增加，這與直接接觸車輪的軟土受壓變形后側(cè)向擠壓相鄰軟土產(chǎn)生的應(yīng)力有關(guān)。

圖7 不同載重時軟土表面壓應(yīng)力Fig.7 Stress of soft-soil surface under different wheel loads

不同載重條件下，車輪侵入軟土最大深度（umax）和最大壓應(yīng)力（σmax）變化曲線如圖8 所示。隨車輪載重增大（1 000 kg～11 000 kg），車輪侵入軟土最大深度分別為1.59 cm、4.61 cm、9.72 cm、16.19 cm、24.01 cm 和31.27 cm。載重增大至11 倍，侵入深度增大至19.67 倍，侵入深度相對改變量約為179%。當(dāng)載重小于5 000 kg 時，車輪侵入軟土最大深度隨載重非線性增加。當(dāng)載重大于5 000 kg 時，車輪侵入軟土最大深度幾乎隨載重線性增加。此外，當(dāng)載重小于5 000 kg時，最大壓應(yīng)力幾乎與載重線性相關(guān)。當(dāng)載重大于5 000 kg 時，最大壓應(yīng)力隨載重非線性增加。隨侵入深度增加，車輪-軟土接觸面積也持續(xù)增加，但最大壓應(yīng)力增加趨勢并未改變，說明載重是侵入深度的關(guān)鍵參數(shù)之一。同時由于侵入深度與載重的關(guān)系是單調(diào)正相關(guān)且非線性不強(qiáng)，因而在現(xiàn)場車輛通行性快速評估過程中載重變化引起的侵入深度比較容易預(yù)測，并可近似按線性擬合。

圖8 最大侵入深度-載重-最大壓應(yīng)力關(guān)系Fig.8 Relationship between maximum intrusion depth，loads，and maximum stress

4.2 軟土彈性模量的影響

改變軟土彈性模量（4 MPa、8 MPa、12 MPa、16 MPa、20 MPa 和24 MPa），分析車輪-軟土相互作用規(guī)律。選取車輪徑向?qū)ΨQ面在軟土表面的投影線及其延伸線為特征線，車輪中心的投影點(diǎn)為0 m 位置。彈性模量變化后，特征線上的軟土表面壓應(yīng)力曲線如圖9 所示。不同彈性模量時，車輪對軟土的壓應(yīng)力均以車輪中心投影點(diǎn)為對稱點(diǎn)接近對稱分布，對稱點(diǎn)（車輪中心的投影點(diǎn)）壓應(yīng)力最大，壓應(yīng)力輪廓線幾乎不變，說明車輪接觸軟土的面積也幾乎不變。隨彈性模量載重增加，壓應(yīng)力最大值逐漸降低，分別為51.80 kPa、48.90 kPa、48.55 kPa、48.40 kPa、48.30 kPa 和48.20 kPa。車輪徑向壓應(yīng)力減增拐點(diǎn)位置幾乎不變，說明與車輪直接接觸的軟土變形有限，相鄰的軟土受側(cè)向擠壓程度相似。

圖9 不同彈性模量時軟土表面壓應(yīng)力Fig.9 Stress of soft-soil surface with different elastic modulus

不同彈性模量條件下，車輪侵入軟土最大深度（umax）和最大壓應(yīng)力（σmax）變化曲線如圖10 所示。隨彈性模量增大（4 MPa～24 MPa），車輪侵入軟土最大深度分別為14.68 cm、13.02 cm、11.92 cm、11.31 cm、10.95 cm 和10.83 cm。彈性模量增至6 倍，侵入深度約降低至71.43%，侵入深度相對改變量約為23%。車輪侵入軟土最大深度和最大壓應(yīng)力均隨彈性模量增大非線性降低。當(dāng)彈性模量大于8 MPa，最大壓應(yīng)力幾乎不變，即車輪-軟土地面接觸面積也幾乎不變，說明軟土處于彈性變形階段，尚未發(fā)生塑性變形。當(dāng)彈性模量較大時（＞10 MPa），增大彈性模量后，最大侵入深度降低約9%。當(dāng)載荷保持一定時，對處于彈性變形階段的軟土，彈性模量變化對車輪侵入深度影響有限。同時由于侵入深度與彈性模量的關(guān)系是單調(diào)負(fù)相關(guān)且變化量不大，因而在現(xiàn)場車輛通行性快速評估過程中彈性模量對侵入深度的影響比較容易預(yù)測。

圖10 最大侵入深度-彈性模量-最大壓應(yīng)力關(guān)系Fig.10 Relationship between maximum intrusion depth，modulus，and maximum stress

4.3 軟土強(qiáng)度的影響

車輪侵入軟土過程中，軟土逐漸從彈性變形階段進(jìn)入塑性變形階段，侵入深度與軟土塑性變形密切相關(guān)，因而軟土強(qiáng)度參數(shù)對侵入深度有較大影響。改變軟土強(qiáng)度參數(shù)，黏聚力c（5 kPa、10 kPa、15 kPa、20 kPa 和25 kPa）和摩擦角φ（10°、15°、20°、25°和30°），分析車輪-軟土相互作用規(guī)律。選取車輪徑向?qū)ΨQ面在軟土表面的投影線及其延伸線為特征線，車輪中心的投影點(diǎn)為0 m位置。以摩擦角10°為例，改變黏聚力（5 kPa、10 kPa、15 kPa、20 kPa 和25 kPa），特征線上的軟土表面壓應(yīng)力曲線如圖11 所示。不同黏聚力時，車輪對軟土的壓應(yīng)力均以車輪中心投影點(diǎn)為對稱點(diǎn)接近對稱分布，對稱點(diǎn)（車輪中心的投影點(diǎn)）壓應(yīng)力最大，壓應(yīng)力輪廓線逐漸縮小，說明車輪接觸軟土的面積逐漸減小。隨黏聚力增加，壓應(yīng)力最大值逐漸增大，分別為48.90 kPa、58.87 kPa、62.39 kPa、63.21 kPa 和63.27 kPa。車輪徑向壓應(yīng)力減增拐點(diǎn)位置與黏聚力大小負(fù)相關(guān)，黏聚力越大拐點(diǎn)距離中心點(diǎn)越近，說明車輪相鄰軟土受側(cè)向擠壓的變形范圍逐漸減小。

圖11 不同黏聚力時軟土表面壓應(yīng)力Fig.11 Stress of soft-soil surface with different cohesion

當(dāng)摩擦角為10°時，改變黏聚力大小，車輪侵入軟土最大深度（umax）和最大壓應(yīng)力（σmax）變化曲線如圖12 所示。不同黏聚力條件下，車輪侵入軟土最大深度分別為13.02 cm、5.83 cm、4.42 cm、3.80 cm 和3.36 cm。黏聚力增至5倍，侵入深度約降低至25.77%，侵入深度相對改變量約為78%。車輪侵入軟土最大深度隨黏聚力增大非線性降低，黏聚力越小，侵入深度對黏聚力越敏感。最大壓應(yīng)力隨黏聚力增大先增加然后趨于穩(wěn)定，即黏聚力改變后車輪-軟土接觸面積逐漸降低至定值。

圖12 最大侵入深度-黏聚力-最大壓應(yīng)力關(guān)系Fig.12 Relationship between maximum intrusion depth，cohesion，and maximum stress

以黏聚力5 kPa 為例，改變摩擦角（10°、15°、20°、25°和30°），特征線上的軟土表面壓應(yīng)力曲線如圖13 所示。不同摩擦角時，車輪對軟土的壓應(yīng)力均以車輪中心投影點(diǎn)為對稱點(diǎn)接近對稱分布，對稱點(diǎn)（車輪中心的投影點(diǎn)）壓應(yīng)力最大。隨摩擦角增加，壓應(yīng)力最大值逐漸增大，分別為48.90 kPa、55.38 kPa、62.79 kPa、69.59 kPa 和69.89 kPa。車輪徑向壓應(yīng)力減增拐點(diǎn)位置與摩擦角大小負(fù)相關(guān)，摩擦角越大拐點(diǎn)距離中心點(diǎn)越近，說明車輪相鄰軟土受側(cè)向擠壓的變形范圍逐漸減小。

圖13 不同摩擦角時軟土表面壓應(yīng)力Fig.13 Stress of soft-soil surface with different friction angles

當(dāng)黏聚力為5 kPa 時，改變摩擦角大小，車輪侵入軟土最大深度（umax）和最大壓應(yīng)力（σmax）變化曲線如圖14 所示。不同摩擦角時，車輪侵入軟土最大深度分別為13.02 cm、6.39 cm、4.15 cm、3.00 cm 和2.40 cm。摩擦角增至3 倍，侵入深度約降至18.42%，侵入深度相對改變量為181%。車輪侵入軟土最大深度隨摩擦角增大非線性降低，摩擦角越小，侵入深度對摩擦角越敏感。最大壓應(yīng)力隨摩擦角增大先增加然后保持一定。

圖14 最大侵入深度-摩擦角-最大壓應(yīng)力關(guān)系曲線Fig.14 Relationship between maximum intrusion depth，friction angle，and maximum stress

當(dāng)黏聚力變化或摩擦角單獨(dú)變化時，最大侵入深度均非線性單調(diào)變化。本文中摩擦角對車輪侵入深度的影響更大，這是由于壓應(yīng)力大于黏聚力，軟土抗剪強(qiáng)度與黏聚力呈線性關(guān)系，而與摩擦角為正切函數(shù)關(guān)系。當(dāng)壓應(yīng)力和黏聚力處于同一量級時，侵入深度主要受黏聚力和摩擦角共同影響，且對摩擦角更為敏感。當(dāng)壓應(yīng)力量級高于黏聚力時，侵入深度主要受摩擦角影響。當(dāng)壓應(yīng)力量級低于黏聚力時，侵入深度主要受黏聚力影響。這種雙因素影響比載重或彈性模量單因素影響更為復(fù)雜。進(jìn)一步，黏聚力和摩擦角共同決定的軟土抗剪強(qiáng)度可能更為直接的影響車輪侵入深度。準(zhǔn)靜態(tài)條件下，分析車輪侵入深度與軟土抗剪強(qiáng)度（τmax=c+σmaxtanφ）的關(guān)系，如圖15 所示。隨軟土抗剪強(qiáng)度的增加，車輪侵入深度非線性降低。當(dāng)軟土最大剪應(yīng)力較低時，車輪侵入深度有顯著變化，這與單獨(dú)考慮黏聚力或摩擦角時的規(guī)律相似。軟土抗剪強(qiáng)度增至3.31 倍，侵入深度約降低至18.45%，侵入深度相對改變量約為164%，與摩擦角的影響處于同一量級。

圖15 侵入深度-最大剪應(yīng)力關(guān)系Fig.15 Relationship between intrusion depth and maximum shear stress

裝備通行性評估中，通常需要快速得到不同軟土地面上車輪的侵入深度。軟土強(qiáng)度特性對侵入深度的非線性影響非常明顯，但抗剪強(qiáng)度往往不容易在現(xiàn)場直接獲得，而黏聚力和摩擦角可以通過現(xiàn)場快速勘測解譯獲得。通過圖12、圖14和圖15 分析可得，侵入深度與抗剪強(qiáng)度的關(guān)系與其單獨(dú)考慮黏聚力或摩擦角的關(guān)系相似。因此，當(dāng)現(xiàn)場軟土抗剪強(qiáng)度不容易獲得時，可以考慮通過黏聚力和摩擦角的變化，間接獲得車輪侵入深度。以車輪載荷為6 000 kg 例，改變黏聚力（5 kPa、10 kPa、15 kPa、20 kPa 和25 kPa）和摩擦角（10°、15°、20°、25°和30°），建立黏聚力、摩擦角和侵入深度構(gòu)成的三維曲面，如圖16 所示。侵入深度最大值和最小值在曲面對角線的兩端，最大值附近區(qū)域侵入深度梯度較大。當(dāng)黏聚力大于10 kPa 且摩擦角大于15°時，侵入深度變化有限。因此，基于載重與侵入深度近似線性的關(guān)系，建立某型裝備不同載重條件下的軟土侵入深度曲面數(shù)據(jù)庫，可為快速確定現(xiàn)場測量軟土關(guān)鍵力學(xué)參數(shù)提供理論依據(jù)，有助于優(yōu)化裝備通行性評估策略，提高裝備通行性評估效率。

圖16 黏聚力-摩擦角-侵入深度關(guān)系曲面Fig.16 Relationship between cohesion，friction angle and intrusion depth

5 結(jié)論

本文提出了CDEM-MPM 耦合計算方法，利用CDEM 模擬車輪，利用MPM 模擬軟土，可有效模擬彈性車輪侵入軟土過程中的大變形問題，同時降低了數(shù)值模擬所需的計算成本。通過建立車輪-軟土地面耦合模型，采用CDEM-MPM耦合計算方法，數(shù)值模擬了車輪和軟土地面相互作用過程，定量分析了車輪載重、軟土彈性模量、軟土強(qiáng)度參數(shù)（黏聚力、摩擦角）與軟土表面壓應(yīng)力、侵入深度等參數(shù)之間的關(guān)系，獲得了車輪侵入軟土深度的變化規(guī)律。結(jié)論如下：

（1）當(dāng)載重增大11 倍，侵入深度增大至19.67 倍，侵入深度相對改變量約為179%。載重是影響侵入深度的關(guān)鍵參數(shù)之一。車輪侵入深度與車輪載重接近線性正相關(guān)，因而在現(xiàn)場車輛通行性快速評估過程中載重引起的侵入深度變化比較容易預(yù)測，并可近似按線性擬合。

（2）當(dāng)彈性模量增至6 倍，侵入深度降低至71.43%，侵入深度相對改變量約為23%。同時，隨彈性模量持續(xù)增加，侵入深度值降低量小于9%。當(dāng)載荷保持一定時，對處于彈性變形階段的軟土，侵入深度值變化不大。車輪侵入深度對彈性模量不敏感，因而在現(xiàn)場車輛通行性快速評估過程中彈性模量引起的侵入深度變化也比較容易預(yù)測。

（3）當(dāng)黏聚力增至5 倍，侵入深度降低至25.77%，侵入深度相對改變量約為78%。摩擦角增至3 倍，侵入深度降低至18.42%，侵入深度相對改變量約為181%。黏聚力和摩擦角共同影響的抗剪強(qiáng)度增至3.31 倍，侵入深度降低至18.5%，侵入深度相對改變量約為164%。侵入深度受到兩個強(qiáng)度參數(shù)的共同影響，且對摩擦角更為敏感，這導(dǎo)致侵入深度較難實現(xiàn)現(xiàn)場預(yù)測。

基于CDEM-MPM 耦合數(shù)值計算方法，建立某型號裝備不同載重條件下的軟土侵入深度三維曲面數(shù)據(jù)庫，可為軟土地面裝備通行性評估中現(xiàn)場待測的關(guān)鍵力學(xué)參數(shù)提供理論依據(jù)。此外，基于該數(shù)值計算方法，下一步還可分析車輪不同運(yùn)動狀態(tài)與軟土侵入深度的關(guān)系，有助于優(yōu)化裝備通行性評估策略，提高裝備通行性評估效率。