基于CDEM-MPM 耦合的軟地面車輛侵入規(guī)律研究

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

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

1 引言

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

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

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

2 數(shù)值方法

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

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

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

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

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

節(jié)點運動計算如下：

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

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

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

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

局部坐標系下，塊體間的接觸力可以表達為：

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

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

2.2 物質點法（MPM）

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

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

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

離散顆粒的質量密度可以表示為：

其中：Mp為質量，np為顆粒總數(shù)，δ為狄拉克函數(shù)（Dirac Delta Function），xpi為顆粒p的坐標。

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

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

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

其中：uip、uiI、uip，j分別為顆粒位移、網(wǎng)格位移和顆粒位移的導數(shù)。g為有限元的網(wǎng)格節(jié)點總數(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分別表示距離和顆粒半徑。

當CDEM 單元與MPM 物質點接觸時，如圖2（b）所示，法向接觸力和剪切接觸力計算分別如下：

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

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

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

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

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

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

3.1 幾何模型

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

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

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

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

選取某黏性土作為典型軟土地面，該黏性土物理力學參數(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°）時，車輪侵入軟土地面的深度。軟土地面力學參數(shù)變化涵蓋了《工程勘察通用規(guī)范》（GB 55017-2021）中黏性土、粉土和部分砂土的力學參數(shù)取值范圍［26］。計算過程中，每次變化一個參數(shù)，其余參數(shù)按表1 中數(shù)值取值。當車輪氣壓達標時，車輪變形量比軟土變形量低二到三個量級，車輪可采用線彈性本構模型。當車輪變形遠小于軟土變形，可簡化車輪中橡膠、鋼絲、簾布以及帶束等材料各自的力學特性，采用加權平均后整體車輪的力學特性，彈性模量大于470 MPa，泊松比約為0.48。

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

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

4 計算結果

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所示。當載重為3 000 kg 時，車輪侵入軟土深度大約為4.61 cm。當載重為11 000 kg 時，車輪侵入軟土深度大約為31.27 cm。隨載重增加，車輪侵入深度逐漸增大，車輪徑向（進行方向）和軸向軟土受擠壓向上隆起。軟土隆起范圍和高度與載重正相關，車輪侵入軟土過程中，三維效應非常明顯。

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

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

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

不同載重條件下，車輪侵入軟土最大深度（umax）和最大壓應力（σ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%。當載重小于5 000 kg 時，車輪侵入軟土最大深度隨載重非線性增加。當載重大于5 000 kg 時，車輪侵入軟土最大深度幾乎隨載重線性增加。此外，當載重小于5 000 kg時，最大壓應力幾乎與載重線性相關。當載重大于5 000 kg 時，最大壓應力隨載重非線性增加。隨侵入深度增加，車輪-軟土接觸面積也持續(xù)增加，但最大壓應力增加趨勢并未改變，說明載重是侵入深度的關鍵參數(shù)之一。同時由于侵入深度與載重的關系是單調(diào)正相關且非線性不強，因而在現(xiàn)場車輛通行性快速評估過程中載重變化引起的侵入深度比較容易預測，并可近似按線性擬合。

圖8 最大侵入深度-載重-最大壓應力關系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ī)律。選取車輪徑向對稱面在軟土表面的投影線及其延伸線為特征線，車輪中心的投影點為0 m 位置。彈性模量變化后，特征線上的軟土表面壓應力曲線如圖9 所示。不同彈性模量時，車輪對軟土的壓應力均以車輪中心投影點為對稱點接近對稱分布，對稱點（車輪中心的投影點）壓應力最大，壓應力輪廓線幾乎不變，說明車輪接觸軟土的面積也幾乎不變。隨彈性模量載重增加，壓應力最大值逐漸降低，分別為51.80 kPa、48.90 kPa、48.55 kPa、48.40 kPa、48.30 kPa 和48.20 kPa。車輪徑向壓應力減增拐點位置幾乎不變，說明與車輪直接接觸的軟土變形有限，相鄰的軟土受側向擠壓程度相似。

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

不同彈性模量條件下，車輪侵入軟土最大深度（umax）和最大壓應力（σ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%。車輪侵入軟土最大深度和最大壓應力均隨彈性模量增大非線性降低。當彈性模量大于8 MPa，最大壓應力幾乎不變，即車輪-軟土地面接觸面積也幾乎不變，說明軟土處于彈性變形階段，尚未發(fā)生塑性變形。當彈性模量較大時（＞10 MPa），增大彈性模量后，最大侵入深度降低約9%。當載荷保持一定時，對處于彈性變形階段的軟土，彈性模量變化對車輪侵入深度影響有限。同時由于侵入深度與彈性模量的關系是單調(diào)負相關且變化量不大，因而在現(xiàn)場車輛通行性快速評估過程中彈性模量對侵入深度的影響比較容易預測。

圖10 最大侵入深度-彈性模量-最大壓應力關系Fig.10 Relationship between maximum intrusion depth，modulus，and maximum stress

4.3 軟土強度的影響

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

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

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

圖12 最大侵入深度-黏聚力-最大壓應力關系Fig.12 Relationship between maximum intrusion depth，cohesion，and maximum stress

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

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

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

圖14 最大侵入深度-摩擦角-最大壓應力關系曲線Fig.14 Relationship between maximum intrusion depth，friction angle，and maximum stress

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

圖15 侵入深度-最大剪應力關系Fig.15 Relationship between intrusion depth and maximum shear stress

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

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

5 結論

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

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

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

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

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