宰德志，龐 銳

(大連理工大學(xué) 水利工程學(xué)院， 遼寧 大連 116024)

邊坡穩(wěn)定問(wèn)題一直是巖土領(lǐng)域多年來(lái)的熱點(diǎn)研究課題之一。數(shù)值模擬是明確滑坡失穩(wěn)過(guò)程，揭示致災(zāi)機(jī)理的有效手段[1]。有限元法作為一種傳統(tǒng)連續(xù)介質(zhì)方法，常用于判斷邊坡的極限狀態(tài)，但在處理失穩(wěn)后的大變形時(shí)會(huì)遇到網(wǎng)格畸變問(wèn)題[2]。離散元法雖然能夠處理巖土大變形問(wèn)題，但其內(nèi)部參數(shù)的確定難度大，且進(jìn)行大規(guī)模計(jì)算耗時(shí)長(zhǎng)、效率低[3]。

Sulsky等[4]改進(jìn)FLIP方法，并將其擴(kuò)展到固體力學(xué)領(lǐng)域，建立了物質(zhì)點(diǎn)法。為解決標(biāo)準(zhǔn)物質(zhì)點(diǎn)法中質(zhì)點(diǎn)穿越背景網(wǎng)格引起的應(yīng)力振蕩等問(wèn)題，Bardenhagen等[5]對(duì)物質(zhì)點(diǎn)法進(jìn)行擴(kuò)展，提出了廣義插值物質(zhì)點(diǎn)法(GIMP)。國(guó)內(nèi)一些學(xué)者對(duì)物質(zhì)點(diǎn)法的發(fā)展也做出突出貢獻(xiàn)。Zhang等[6]把物質(zhì)點(diǎn)法與有限元法耦合在一起，提出了物質(zhì)點(diǎn)有限元法。Wang等[7]提出了隱式物質(zhì)點(diǎn)法，改進(jìn)算法精度，并將其應(yīng)用在巖土領(lǐng)域中。

物質(zhì)點(diǎn)法是一種完全的拉格朗日質(zhì)點(diǎn)類方法，該方法充分發(fā)揮拉格朗日法和歐拉法各自的優(yōu)點(diǎn)，能夠準(zhǔn)確、高效地模擬巖土材料的大變形行為。本文采用物質(zhì)點(diǎn)開(kāi)展邊坡失穩(wěn)狀態(tài)和滑坡的全過(guò)程模擬，并詳細(xì)分析土體強(qiáng)度對(duì)滑動(dòng)距離的影響，為邊坡災(zāi)害的防治提供參考。

1 物質(zhì)點(diǎn)法基本理論

物質(zhì)點(diǎn)法采用拉格朗日質(zhì)點(diǎn)和歐拉網(wǎng)格來(lái)雙重描述離散域，將連續(xù)體離散成攜帶材料信息的一系列質(zhì)點(diǎn)，如圖1所示。在每一計(jì)算時(shí)步內(nèi)，質(zhì)點(diǎn)和計(jì)算網(wǎng)格固連，避免了數(shù)值計(jì)算中非線性對(duì)流項(xiàng)的產(chǎn)生；計(jì)算時(shí)步結(jié)束后，丟棄已變形的計(jì)算網(wǎng)格，采用新的未變形計(jì)算網(wǎng)格，從而避免網(wǎng)格畸變。其中網(wǎng)格僅用于動(dòng)量方程的求解和空間導(dǎo)數(shù)的計(jì)算，質(zhì)點(diǎn)攜帶該區(qū)域網(wǎng)格的所有信息[8]。

圖1 物質(zhì)點(diǎn)法示意圖

1.1 運(yùn)動(dòng)方程

(1) 動(dòng)量方程[9]：

(1)

(2) 引入虛位移，并在連續(xù)體域內(nèi)積分，得到更新拉格朗日格式的弱形式[8]：

(2)

(3) 將連續(xù)體密度近似，把積分弱形式轉(zhuǎn)化為求和的離散形式[8]：

(3)

式中:np為物質(zhì)點(diǎn)總數(shù)；mp為物質(zhì)點(diǎn)p的質(zhì)量；h為引入的假想邊界層厚度；帶有下標(biāo)p的物理量表示物質(zhì)點(diǎn)p的變量。

(4) 依據(jù)牛頓第二定律，分類、結(jié)合、化簡(jiǎn)得到背景網(wǎng)格結(jié)點(diǎn)的運(yùn)動(dòng)方程[8]：

(4)

1.2 方程求解

標(biāo)準(zhǔn)物質(zhì)點(diǎn)法的求解格式有三種[8]，分別為USF、USL和MUSL。根據(jù)應(yīng)力更新時(shí)刻的差異，劃分為USF和USL兩種求解格式，前者在每個(gè)計(jì)算時(shí)步開(kāi)始時(shí)更新應(yīng)力，后者在計(jì)算時(shí)步結(jié)束時(shí)更新應(yīng)力。MUSL是基于USL改進(jìn)的一種求解格式，將更新的質(zhì)點(diǎn)動(dòng)量映射回背景網(wǎng)格后再計(jì)算節(jié)點(diǎn)速度。MUSL具有良好的能量守恒性，故本文采用MUSL求解格式。

當(dāng)網(wǎng)格尺寸較小時(shí)，標(biāo)準(zhǔn)物質(zhì)點(diǎn)法中質(zhì)點(diǎn)會(huì)穿越背景網(wǎng)格，引起數(shù)值噪音，產(chǎn)生應(yīng)力振蕩[10]。相比于標(biāo)準(zhǔn)物質(zhì)點(diǎn)法，廣義插值物質(zhì)點(diǎn)法能有效減輕應(yīng)力振蕩，但會(huì)顯著增加計(jì)算時(shí)間，且隨著質(zhì)點(diǎn)數(shù)目的增加，計(jì)算時(shí)間增加的趨勢(shì)變快?；谟?jì)算效率考慮，本文的計(jì)算方法采用標(biāo)準(zhǔn)物質(zhì)點(diǎn)法。針對(duì)運(yùn)動(dòng)方程是關(guān)于時(shí)間的二階常微分方程，本文利用中心差分法[11]求解，計(jì)算時(shí)間步長(zhǎng)小于臨界時(shí)間步長(zhǎng)，具體要求滿足式(5)[8]。

(5)

式中：α為時(shí)間步長(zhǎng)因子，本文取0.8，Δl為背景網(wǎng)格長(zhǎng)度；cp為壓縮波波速，與巖土材料特性有關(guān)。

1.3 本構(gòu)模型

(6)

Mohr-Coulomb(MC)屈服面在π平面上為不規(guī)則的六邊形，在六個(gè)頂點(diǎn)處法線方向不唯一，導(dǎo)致頂點(diǎn)處的數(shù)值求解困難。本文采用的Drucker-Prager(DP)屈服面在主應(yīng)力空間中為圓錐面，在π平面上為圓形，便于求解，其屈服函數(shù)滿足式(7)[12]：

(7)

式中：J2為第二偏應(yīng)力張量不變量；I1為第一應(yīng)力張量不變量；qφ與kφ為材料常數(shù)。本文采取DP屈服面在π平面上內(nèi)接MC屈服面，見(jiàn)圖2。材料常數(shù)與黏聚力c及內(nèi)摩擦角φ的關(guān)系如式(8)[12]。

圖2 π平面上MC屈服面和DP屈服面的關(guān)系圖[8]

(8)

在DP模型中，剪切屈服采用非關(guān)聯(lián)塑性流動(dòng)法則，其勢(shì)函數(shù)Ψs的關(guān)系滿足式(9)[8]：

(9)

式中：qΨ為材料常數(shù)，由剪脹角ψ確定，二者的關(guān)系和qφ與φ之間的關(guān)系相同。由于材料滿足塑性不可壓縮條件，剪脹角取為0。拉伸屈服采用關(guān)聯(lián)流動(dòng)法則，其勢(shì)函數(shù)Ψt的表達(dá)式滿足式(10)[8]：

(10)

2 程序驗(yàn)證

Bui等[13]通過(guò)堆積干燥鋁棒坍塌試驗(yàn)?zāi)M砂土的失穩(wěn)過(guò)程，并采用光滑粒子流體動(dòng)力學(xué)(SPH)方法對(duì)其進(jìn)行數(shù)值模擬，同時(shí)對(duì)黏土失穩(wěn)的大變形過(guò)程開(kāi)展了分析。本節(jié)應(yīng)用物質(zhì)點(diǎn)法模擬上述試驗(yàn)和粘土的大變形過(guò)程。程序采用清華大學(xué)張雄教授課題組開(kāi)發(fā)的MPM3D-F90[8]，并在其基礎(chǔ)上進(jìn)行了改編，編譯環(huán)境為Microsoft Visual Studio 2010，數(shù)據(jù)后處理采用Tecplot 360 EX 2018 R1。

2.1 干燥鋁棒坍塌試驗(yàn)?zāi)M驗(yàn)證

為了便于對(duì)比試驗(yàn)結(jié)果，模型尺寸及材料參數(shù)和Bui等[13]的數(shù)據(jù)保持一致，即模型長(zhǎng)200 mm，寬2 mm，高100 mm，材料參數(shù)如表1所示。模型采用40 000個(gè)物質(zhì)點(diǎn)離散，剛體邊界采用10 128個(gè)物質(zhì)點(diǎn)離散，模型和剛體邊界之間采用剛?cè)峤佑|。相鄰物質(zhì)點(diǎn)間距為1 mm，背景網(wǎng)格尺寸為2 mm。模型左(X=-8 mm)、底(Y=-8 mm)邊界采用固定邊界，右(X=500 mm)、上(Y=120 mm)邊界采用自由邊界，平面法向兩邊界(Z=0 mm、Z=2 mm)為對(duì)稱邊界。

表1 干燥鋁棒材料參數(shù)

圖3(a)為Bui等通過(guò)試驗(yàn)得到的干燥鋁棒坍塌最終構(gòu)型圖，圖3(b)為利用物質(zhì)點(diǎn)法數(shù)值模擬得到的最終堆積形態(tài)。可以看出，利用物質(zhì)點(diǎn)法模擬的堆積物最終形態(tài)及滑動(dòng)區(qū)域與非滑動(dòng)區(qū)域的分界線和試驗(yàn)結(jié)果均吻合較好，從而驗(yàn)證了該程序模擬砂土大變形問(wèn)題的有效性。

2.2 黏土大變形模擬驗(yàn)證

黏土模型長(zhǎng)4 m，寬0.02 m，高2 m。土體材料參數(shù)與Bui等[13]數(shù)值模擬中保持一致，具體參數(shù)如表2所示。模型采用40 000個(gè)物質(zhì)點(diǎn)離散，剛體邊界采用5 728個(gè)物質(zhì)點(diǎn)離散，相鄰物質(zhì)點(diǎn)間距為0.02 m，背景網(wǎng)格尺寸為0.04 m。模型左(X=-0.16 m)、底(Y=-0.16 m)邊界采用固定邊界，平面法向兩邊界(Z=0 m、Z=0.02 m)采用對(duì)稱邊界，右(X=5 m)、上(Y=2 m)邊界為自由邊界。

圖3 鋁棒坍塌堆積形態(tài)圖

表2 黏土材料參數(shù)

圖4(a)圖為Bui等采用SPH方法進(jìn)行數(shù)值模擬得到的黏土塑性偏應(yīng)變?cè)茍D，圖4(b)圖為采用物質(zhì)點(diǎn)法進(jìn)行數(shù)值模擬得到的累計(jì)等效塑性應(yīng)變?cè)茍D。從圖中可得，滑動(dòng)帶的位置、形狀、寬度及穩(wěn)定后的最終構(gòu)型基本保持一致。同時(shí)，兩種方法模擬黏土失穩(wěn)的過(guò)程也具有較高相似性，即相同計(jì)算時(shí)刻，黏土失穩(wěn)形態(tài)相似。這表明該程序能夠合理地模擬黏土大變形問(wèn)題。

圖4 SPH模擬圖[13]與MPM模擬圖

3 影響因素分析

邊坡滑動(dòng)距離是指坡腳到失穩(wěn)后滑坡前緣的距離，是邊坡安全控制的重要因素之一。針對(duì)邊坡的滑動(dòng)失穩(wěn)，分析了不同坡角(30°、45°、60°、75°)情況下邊坡的內(nèi)摩擦角和黏聚力對(duì)滑動(dòng)距離的影響規(guī)律。邊坡幾何尺寸如圖5(以45°坡角為例)所示，底部為固定邊界，上部為自由邊界，其余為對(duì)稱邊界。坡體材料參數(shù)匯總于表3[8]。

圖5 邊坡幾何尺寸

表3 邊坡土體參數(shù)

基于滑動(dòng)距離和邊坡形態(tài)，以45°坡角的邊坡為例，確定合理的相鄰物質(zhì)點(diǎn)間距[14]。邊坡模型見(jiàn)圖5，內(nèi)摩擦角取20°，黏聚力為5.0 kPa。采用2 460、9 820、15 475和39 240個(gè)物質(zhì)點(diǎn)離散連續(xù)體，其對(duì)應(yīng)的相鄰物質(zhì)點(diǎn)間距分別為1.00 m、0.50 m、0.40 m和0.25 m。數(shù)值模擬得到的構(gòu)型圖見(jiàn)圖6，圖中顏色差異代表不同的等效塑性應(yīng)變值。采用上述物質(zhì)點(diǎn)數(shù)進(jìn)行數(shù)值模擬得到的滑動(dòng)距離分別為8.0 m、16.6 m、17.4 m、18.1 m，計(jì)算耗時(shí)分別為83 s、340 s、678 s、6 848 s(基于2.9GHzCPU(Intel(R) Core(TM)i7-10700 CPU)和16 GB內(nèi)存計(jì)算機(jī))。考慮邊坡最終構(gòu)型、滑動(dòng)距離及計(jì)算耗時(shí)，采用15 475個(gè)物質(zhì)點(diǎn)進(jìn)行數(shù)值模擬，能夠保證結(jié)果的精度且計(jì)算效率高。因此，在后續(xù)參數(shù)研究中，相鄰物質(zhì)點(diǎn)間距取0.4 m。

3.1 內(nèi)摩擦角

固定黏聚力值為5 kPa，在10°～35°范圍內(nèi)取不同內(nèi)摩擦角，分析不同坡角情況下，內(nèi)摩擦角與滑動(dòng)距離的關(guān)系，結(jié)果如圖7所示?？梢钥闯觯吰率Х€(wěn)后其滑動(dòng)距離與土體內(nèi)摩擦角呈現(xiàn)負(fù)相關(guān)性。坡角為30°的邊坡在內(nèi)摩擦角增至25°時(shí)滑動(dòng)距離為0，而坡角為45°的邊坡在內(nèi)摩擦角增至35°時(shí)滑動(dòng)距離為0。內(nèi)摩擦角較小(<25°)時(shí)，滑動(dòng)距離對(duì)其變化比較敏感；內(nèi)摩擦角較大時(shí)，滑動(dòng)距離的變化程度減小，趨于常值。內(nèi)摩擦角一定時(shí)，邊坡越陡，失穩(wěn)后的滑動(dòng)距離越大，其可能造成的危害愈大。數(shù)值模擬過(guò)程中，內(nèi)摩擦角10°時(shí)，滑裂面底部低于坡腳高度，失穩(wěn)模式類似于深層滑坡，且坡腳以上邊坡發(fā)生巨大位移；內(nèi)摩擦角35°時(shí)，陡坡滑動(dòng)距離較小，緩坡無(wú)滑動(dòng)，只是在邊坡的中上部發(fā)生較大塑性應(yīng)變。同時(shí)，失穩(wěn)邊坡達(dá)到穩(wěn)定狀態(tài)后，邊坡傾角隨內(nèi)摩擦角增大逐漸增加，但始終小于后者。

圖6 不同物質(zhì)點(diǎn)數(shù)目邊坡形態(tài)圖

圖7 內(nèi)摩擦角與滑動(dòng)距離關(guān)系圖

3.2 黏聚力

固定內(nèi)摩擦角值為20°，在0 kPa～24 kPa范圍內(nèi)取不同黏聚力，分析不同坡角情況下，黏聚力與滑動(dòng)距離的關(guān)系，結(jié)果如圖8所示。可以看出，邊坡失穩(wěn)后其滑動(dòng)距離與土體黏聚力呈現(xiàn)負(fù)相關(guān)性。坡角為30°的邊坡在黏聚力增至8 kPa時(shí)滑動(dòng)距離為0，而坡角為45°的邊坡在黏聚力增至24 kPa時(shí)滑動(dòng)距離為0。在黏聚力較小的初始階段(<6 kPa)，黏聚力對(duì)滑動(dòng)距離的影響比較明顯；在黏聚力較大的末尾階段(>18 kPa)，黏聚力對(duì)滑動(dòng)距離的影響趨于穩(wěn)定。黏聚力一定時(shí)，坡角對(duì)滑動(dòng)距離的影響頗大，坡度越大，滑動(dòng)距離越大。數(shù)值模擬過(guò)程中，黏聚力增至24 kPa，形成了坡腳到坡頂?shù)呢灤┗褞В刃苄詰?yīng)變值總體偏小，僅在陡坡中出現(xiàn)較小的滑動(dòng)距離，未發(fā)生明顯的塑性流動(dòng)現(xiàn)象。

圖8 黏聚力與滑動(dòng)距離關(guān)系圖

4 結(jié) 論

(1) 合理地選取離散點(diǎn)數(shù)和間距，對(duì)物質(zhì)點(diǎn)法數(shù)值模擬的結(jié)果至關(guān)重要。建議根據(jù)邊坡最終構(gòu)型、滑動(dòng)距離及計(jì)算效率選取合適的相鄰物質(zhì)點(diǎn)間距。

(2) 滑動(dòng)距離隨土體強(qiáng)度的增大逐漸減小。

(3) 邊坡在地震荷載作用下的穩(wěn)定性也是工程設(shè)計(jì)關(guān)注的重點(diǎn)問(wèn)題，后續(xù)工作將基于物質(zhì)點(diǎn)法深入開(kāi)展邊坡動(dòng)力失穩(wěn)機(jī)理和滑動(dòng)過(guò)程研究。