基于光滑粒子流體力學法分析坡高和坡角對邊坡破壞特征的影響

冉光炯* 劉 勇 白 皓

（四川高速公路建設開發集團有限公司，四川 成都610000）

在我國地形陡峭的山區，由于地質環境極為特殊，地質作用非常復雜，極易導致邊坡失穩破壞[1]。

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

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

1 光滑粒子流體力學法

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

以下簡要介紹土體運動方程的SPH 離散化。連續體的運動可以用以下方程來描述：

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

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

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

公式(4)在有效應力張量已知的情況下，可以使用標準的蛙跳算法[14]進行積分。

本構模型描述了給定材料的應力和應變之間的關系。根據經典塑性理論，將彈塑性材料的總應變率張量ε˙分解為兩部分：彈性應變率張量ε˙e和塑性應變率張量ε˙p：

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

塑性應變率張量ε˙αβp按塑性流動法則計算：

其中λ˙是塑性比例系數的變化率，gp是塑性勢函數。本文采用非關聯塑性流動規律的Drucker-Prager 模型[9]，假設屈服面固定在應力空間中，只有當應力狀態達到屈服面時才會發生塑性變形。因此，只有在滿足下列屈服準則的情況下，才會發生塑性變形：

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

上述土體本構模型需要6 個土體參數：粘聚力系數、摩擦角、剪脹角、楊氏模量E、泊松比和土體密度ρ。Bui 等人[12]在SPH 框架內對當前土體本構模型進行了驗證。其中SPH 關于土體材料重力流動和粘性土上的地基問題的計算結果與實驗，和有限元法一致。

2 模型構建與計算流程

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

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

圖2 Dual-SPH 構建的邊坡模型

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

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

確定了滑動粒子后，在水平方向上用滑動距離S 來量化評價邊坡失穩后的破壞特征，滑動距離S 是指從邊坡破壞前的坡腳到沉積穩定后的邊界點的距離。

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

表2 計算中用到的土體參數取值

表3 25 個邊坡的坡腳和坡高及其滑動距離的計算結果

3 坡高與坡腳對邊坡穩定性的影響分析

本小節在第二小節的基礎上，主要針對邊坡的坡腳和坡高對其穩定性的影響展開分析。首先，采用第2 小節所述的方法，在中分別構建了25 個坡腳和坡高不同的邊坡，具體每個邊坡的坡高和坡腳如表3 所示。然后采用Dual-SPH 分別計算每個邊坡的滑動距離，其結果表3 所示。

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

圖4 邊坡高度、邊坡坡度與滑動距離的關系

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

4 結論

本文基于開源的Dual-SPH 軟件來對邊坡坡高和坡腳對邊坡失穩后滑動距離的影響進行研究，分別構建了25 個不同坡高和坡腳的邊坡，并對其失穩后的滑動距離進行了計算。計算結果表明：隨著邊坡坡度的增大和邊坡高度的增高，邊坡的滑動距離隨之明顯增大。進一步分析發現，在本文所建立的模型與參數條件下，邊坡坡度、坡高與滑動距離的關系可以用二次曲面很好的擬合，且擬合優度R2為99.7%。本文的研究結果為定量分析邊坡坡高、坡腳與邊坡失穩后滑動距離的關系奠定了基礎。后續將繼續針對邊坡參數對其失穩后破壞特征參數的影響展開定量研究。