冉光炯* 劉 勇 白 皓

（四川高速公路建設(shè)開發(fā)集團(tuán)有限公司，四川 成都610000）

在我國地形陡峭的山區(qū)，由于地質(zhì)環(huán)境極為特殊，地質(zhì)作用非常復(fù)雜，極易導(dǎo)致邊坡失穩(wěn)破壞[1]。

在已有的文獻(xiàn)中，研究人員常采用FEM（有限元方法）[2]、DEM（離散元法）[3]和FDM（有限差分法）[4]來研究失穩(wěn)破壞后邊坡的變形特征。然而，在模擬邊坡失穩(wěn)破壞時(shí)，有限元法和有限差分法都具有一定的缺點(diǎn)，如網(wǎng)格變形過大土體破壞引起的變形。DEM法雖然有利于模擬大變形問題，但DEM在處理工程尺度的模擬問題時(shí)，其計(jì)算效率還有待提高，且DEM的微觀參數(shù)難以標(biāo)定。

近年來，無網(wǎng)格平滑粒子流體力學(xué)（SPH）[5-8]已應(yīng)用于探究巖土材料的大變形問題。已有研究表明，該方法可以有效地預(yù)測邊坡破壞引起的土體大變形。本文采用SPH 方法對(duì)該問題進(jìn)行了研究，采用彈塑性Drucker-Prager 模型[9]來模擬土體，探究了不同坡高和坡角的邊坡的破壞特征。

1 光滑粒子流體力學(xué)法

在SPH 方法中，使用一組運(yùn)動(dòng)粒子來模擬連續(xù)介質(zhì)的運(yùn)動(dòng)；每個(gè)都分配有恒定的質(zhì)量，并在相應(yīng)位置“攜帶”場變量。通過加權(quán)求和，將連續(xù)場從各種粒子中插值出來，其中權(quán)重根據(jù)假定的核函數(shù)隨距離而下降。然后，“粒子間相互作用”的范圍對(duì)應(yīng)于內(nèi)核的“支持域”的半徑。類似地，空間導(dǎo)數(shù)是通過在支撐域上的插值來計(jì)算的。SPH 插值理論可參考Liu & Liu[10]和Monaghan[11]的論文。

以下簡要介紹土體運(yùn)動(dòng)方程的SPH 離散化。連續(xù)體的運(yùn)動(dòng)可以用以下方程來描述：

在SPH 框架中，此等式通常離散為

其中i 表示考慮中的粒子；ρi和ρj分別為粒子i 和j 的密度；N 是“相鄰粒子”的數(shù)目，即在粒子i 的支撐域中的粒子的數(shù)目；mj是粒子j 的質(zhì)量；Cαβij是一個(gè)穩(wěn)定項(xiàng)，用于消除應(yīng)力波動(dòng)和拉伸不穩(wěn)定[12]；Wij是核函數(shù)，這里選擇的是三次樣條函數(shù)[13]。公式(4)中的穩(wěn)定項(xiàng)由兩個(gè)分量組成：人工粘度和人工應(yīng)力，其計(jì)算方法與文獻(xiàn)[12]相似。但人工粘性項(xiàng)的聲速在此由下式計(jì)算

式中，E 是土的楊氏模量。

公式(4)在有效應(yīng)力張量已知的情況下，可以使用標(biāo)準(zhǔn)的蛙跳算法[14]進(jìn)行積分。

本構(gòu)模型描述了給定材料的應(yīng)力和應(yīng)變之間的關(guān)系。根據(jù)經(jīng)典塑性理論，將彈塑性材料的總應(yīng)變率張量ε˙分解為兩部分：彈性應(yīng)變率張量ε˙e和塑性應(yīng)變率張量ε˙p：

式中，S˙'αβ為偏有效剪應(yīng)力張量；G 為剪切模量，σ˙'m為平均有效應(yīng)力。

塑性應(yīng)變率張量ε˙αβp按塑性流動(dòng)法則計(jì)算：

其中λ˙是塑性比例系數(shù)的變化率，gp是塑性勢函數(shù)。本文采用非關(guān)聯(lián)塑性流動(dòng)規(guī)律的Drucker-Prager 模型[9]，假設(shè)屈服面固定在應(yīng)力空間中，只有當(dāng)應(yīng)力狀態(tài)達(dá)到屈服面時(shí)才會(huì)發(fā)生塑性變形。因此，只有在滿足下列屈服準(zhǔn)則的情況下，才會(huì)發(fā)生塑性變形：

其中I1和J2是應(yīng)力張量的第一和第二不變量；α?和kc是Drucker-Prager 常數(shù)，這是根據(jù)庫侖材料常數(shù)c(內(nèi)聚力)和φ(內(nèi)摩擦角)計(jì)算得出的。

上述土體本構(gòu)模型需要6 個(gè)土體參數(shù)：粘聚力系數(shù)、摩擦角、剪脹角、楊氏模量E、泊松比和土體密度ρ。Bui 等人[12]在SPH 框架內(nèi)對(duì)當(dāng)前土體本構(gòu)模型進(jìn)行了驗(yàn)證。其中SPH 關(guān)于土體材料重力流動(dòng)和粘性土上的地基問題的計(jì)算結(jié)果與實(shí)驗(yàn)，和有限元法一致。

2 模型構(gòu)建與計(jì)算流程

本文基于開源的Dual-SPH 軟件[15]，采用編寫命令流的方式來進(jìn)行不同坡高和坡腳的邊坡模型的構(gòu)建。其中，坡高和坡腳的定義如圖1 所示，采用Dual-SPH 構(gòu)建出來的邊坡模型如圖2所示。注意在建立模型的時(shí)候，我們將模型沿垂直直面方向拉伸了一個(gè)單位長度(1mm)的寬度，使其形成一個(gè)假三維模型。

圖1 邊坡坡高y(m)與坡腳θ(°)的示意圖

圖2 Dual-SPH 構(gòu)建的邊坡模型

為了探究土體參數(shù)對(duì)邊坡破壞特征的影響規(guī)律，采用滑動(dòng)距離S 來定量表征邊坡的破壞特征。圖2 是通過Dual-SPH 模擬得到的邊坡土體破壞后Dual-SPH 粒子的位移云圖。圖3 用不同的顏色表示不同程度的土體的位移大小，可根據(jù)位移大小變化識(shí)別失穩(wěn)滑動(dòng)帶。

圖3 Dual-SPH 邊坡模型中粒子的位移圖

確定了滑動(dòng)粒子后，在水平方向上用滑動(dòng)距離S 來量化評(píng)價(jià)邊坡失穩(wěn)后的破壞特征，滑動(dòng)距離S 是指從邊坡破壞前的坡腳到沉積穩(wěn)定后的邊界點(diǎn)的距離。

為了控制變量，在所有邊坡模型中，其土體參數(shù)的取值都一致，見表2。

表2 計(jì)算中用到的土體參數(shù)取值

表3 25 個(gè)邊坡的坡腳和坡高及其滑動(dòng)距離的計(jì)算結(jié)果

3 坡高與坡腳對(duì)邊坡穩(wěn)定性的影響分析

本小節(jié)在第二小節(jié)的基礎(chǔ)上，主要針對(duì)邊坡的坡腳和坡高對(duì)其穩(wěn)定性的影響展開分析。首先，采用第2 小節(jié)所述的方法，在中分別構(gòu)建了25 個(gè)坡腳和坡高不同的邊坡，具體每個(gè)邊坡的坡高和坡腳如表3 所示。然后采用Dual-SPH 分別計(jì)算每個(gè)邊坡的滑動(dòng)距離，其結(jié)果表3 所示。

為了便于直觀地觀察邊坡高度與邊坡坡度對(duì)滑動(dòng)距離的影響，我們將邊坡高度、邊坡坡度與滑動(dòng)距離繪制在同一個(gè)三維圖上，如圖4 所示?？梢钥吹剑S著邊坡坡度的增大和邊坡高度的增高，邊坡的滑動(dòng)距離隨之顯著增大。為了研究三者之間的定量關(guān)系，我們用二次曲面對(duì)這些散點(diǎn)進(jìn)行擬合，得到邊坡坡度、邊坡高度與滑動(dòng)距離的關(guān)系式如公式（17）所示，其中x 表示邊坡坡度，y 表示邊坡高度，S 表示邊坡的滑動(dòng)距離。從圖中可以看出，該公式的擬合優(yōu)度R2為99.7%。這表明：在本文的分析框架與參數(shù)條件內(nèi)，邊坡高度、邊坡坡度與滑動(dòng)距離的關(guān)系可以很好地用二次曲面來進(jìn)行擬合。圖5 分別展示了不同邊坡高度與邊坡坡度的SPH 粒子位移云圖。

圖4 邊坡高度、邊坡坡度與滑動(dòng)距離的關(guān)系

圖5 不同邊坡高度與邊坡坡度的SPH 邊坡的位移云圖

4 結(jié)論

本文基于開源的Dual-SPH 軟件來對(duì)邊坡坡高和坡腳對(duì)邊坡失穩(wěn)后滑動(dòng)距離的影響進(jìn)行研究，分別構(gòu)建了25 個(gè)不同坡高和坡腳的邊坡，并對(duì)其失穩(wěn)后的滑動(dòng)距離進(jìn)行了計(jì)算。計(jì)算結(jié)果表明：隨著邊坡坡度的增大和邊坡高度的增高，邊坡的滑動(dòng)距離隨之明顯增大。進(jìn)一步分析發(fā)現(xiàn)，在本文所建立的模型與參數(shù)條件下，邊坡坡度、坡高與滑動(dòng)距離的關(guān)系可以用二次曲面很好的擬合，且擬合優(yōu)度R2為99.7%。本文的研究結(jié)果為定量分析邊坡坡高、坡腳與邊坡失穩(wěn)后滑動(dòng)距離的關(guān)系奠定了基礎(chǔ)。后續(xù)將繼續(xù)針對(duì)邊坡參數(shù)對(duì)其失穩(wěn)后破壞特征參數(shù)的影響展開定量研究。