基于SPH法的排土場邊坡破壞過程數值模擬研究

魏良針，王建平，楊 磊，應濤濤，龔禮岳*，吳則祥

（1.溫州市鐵路與軌道交通投資集團有限公司，浙江 溫州 325000；2.溫州大學 建筑與土木工程學院，浙江 溫州 325035）

邊坡失穩破壞導致的災害具有發生頻率高、突發性強、經濟損失大、社會影響大等特點，排土場邊坡的破壞作為其中典型問題頻發，使得排土場的滑坡災害的防治一直以來都受到較多的關注，因而排土場邊坡的破壞特征分析是邊坡工程中的一個重要研究課題，也一直是巖土工程領域的熱點研究方向之一[1]。

目前對排土場邊坡穩定性的研究主要集中在邊坡失穩的起始階段[2-4]，然而滑坡的整個過程不僅僅只包括滑坡的起始，還包括滑動土體的遷移和最終的堆積。另外，值得注意的是，滑坡的最終后果往往取決于滑動土體的搬運和最終沉積。這可能是由于以往的常規方法在模擬滑坡發生后土體的大變形方面存在困難所致。然而，隨著計算機技術的不斷發展與硬件性能的逐漸提高，近年來發展起來的拉格朗日型無網格粒子法卻能夠很好地解決這類問題，特別適合模擬巖土工程中的大變形問題，例如排土場整體破壞的全過程模擬，而且具有較高的計算效率和較好的數值穩定性。

近年來，基于拉格朗日粒子的無網格方法得到了發展，并越發流行。這些方法在網格法難以建模的問題上更具優勢，因為它們或多或少地避免了網格的使用。常見的方法包括：光滑粒子流體力學法（SPH）[5]，物質點法（MPM）[6]，粒子有限元法（PFEM）[7]。在這些方法中，SPH確立已久，并且在天體物理、流體動力學和固體力學中有著悠久的應用歷史。近年來，SPH還被用于巖土建模領域中。已知的應用包括大型形變，顆粒流動，土地結構間相互作用，滑坡和泥石流、流體-土壤混合物動力學，開裂和地下爆炸[8-11]等一系列方向。結果表明，SPH能夠很好地解決基于網格的方法所不能解決的上述問題，而排土場邊坡滑坡問題的模擬也取得了較大的進步。因此，SPH是有意義的，并有潛力應用于更廣泛的學術和工程巖土應用[12-14]。

在以往學者的研究中，鮮有針對排土場邊坡的破壞特征影響的研究。因此，本文基于SPH（Smoothed particle hydrodynamics）方法，選了25組不同的粘聚力與摩擦角的組合研究其對均質土體邊坡的滑動距離的影響，并在此基礎上模擬了深圳一處排土場的滑坡現象。在邊界條件復雜，破壞面位置和形式未知的條件下，本文基于SPH法模擬了排土場邊坡失穩的完整過程，準確清晰地反映了土坡的破壞機理，計算得出土體內部各處的應力與應變，為此類滑坡災害的防范與治理提供了數據支持。

1 SPH法

SPH是一種基于拉格朗日粒子法，因此，在巖土問題中，SPH控制方程寫成如下拉格朗日形式：

其中，ρ表示密度，ν表示速度，σ表示柯西壓力張量，g表示由外力引起的加速度，如絕大多數巖土問題中常見的重力加速度，表示材料時間導數，? 表示梯度算子。我們沒有考慮孔隙壓力的情況，所以材料可能是干燥的或者可以被總壓力建模。在許多仿真模擬中，方程（1），即連續方程可以刪除，因為我們應用了拉格朗日粒子，質量守恒總是滿足的。連續方程只有在基于狀態方程（EOS）的本構模型中使用。

在SPH中，計算域使用粒子進行劃分，粒子攜帶場變量并隨材料移動。控制方程可以通過跟蹤粒子的運動和攜帶變量的變化來求解。要做到這一點，必須使用SPH插值離散化控制方程。場函數可以由以下積分插值近似：

其中，＜＞表示值的近似，Ω是核函數W（x-x′，h）的影響域，在本文之后簡寫為W。W取決于距離||x-x′||和參數h，即平滑步長。本文的球形支持域W半徑為2h。此外，核函數W必須滿足一定的條件，如歸一化條件，增量函數性質和緊湊支持條件。常用的核函數包括高斯核函數（Gaussian kernel），三次樣條核函數（cubic spine kernel）和溫德蘭核函數（Wendland kernel）。在本文中，由于可以有效阻止粒子聚集，我們只使用C2溫德蘭核函數。

類似的，場函數的空間導數可以通過以下方程近似：

其中，?x′表示導數在x′處取值。上述方程僅在支撐域與計算域邊界不相交時有效。然而，這個條件在邊界附近的區域通常是不滿足的。因此，SPH需要一些邊界的特殊處理來解決所謂的邊界缺陷。

由于可以在支撐域中累加所有粒子的作用，SPH中方程（3）（4）的連續積分插值形式可以重寫為以下形式：

其中，方便起見我們省略了＜＞符號；xi是場函數估計的粒子，xj是支撐域中的粒子；?iWij是 ?xiW（xi-xj，h）的簡寫形式。由于核函數的對稱性，?iWij=- ?jWij；mj表示粒子j的質量，表示粒子的體積。

考慮SPH插值方程（5）（6）的離散形式，控制方程可以重寫為以下形式：

盡管其他文獻中也有控制方程的不同離散形式，上述方程能夠保留線性和角動量，并被廣泛使用。

2 模型構建與邊坡破壞特征量化方法

本文基于C++編寫了SPH計算模擬程序，采用編寫命令流的方式來構建邊坡模型，模型尺寸如圖1所示。所對應的SPH粒子模型如圖2所示。

圖1 邊坡尺寸示意圖

圖2 基于SPH粒子的邊坡模型

對于如圖2所示的二維SPH邊坡，其最靠底面與側面的粒子為受約束粒子，這些粒子的x、y兩個方向的位移都始終為0。

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

圖3 邊坡破壞的SPH粒子位移云圖

確定了滑動粒子后，在水平方向上用Run-out Distance定量描述邊坡的破壞特征，Run-out Distance是指從邊坡破壞前的坡腳到沉積穩定后的邊界點的距離。

3 粘聚力與摩擦角對邊坡破壞特征的影響分析

本節主要針對粘聚力和摩擦角對邊坡破壞特征的影響展開分析。基于第一節介紹的SPH法與第二節所構建的邊坡模型，給邊坡賦值不同的粘聚力和摩擦角進行分析。共對25組不同粘聚力和摩擦角的組合進行了分析，每組粘聚力和摩擦角的取值見表1。其余土體參數的取值均相同，即重度為γ=20 kN/m3，彈性模量為E=20 MPa，泊松比為μ=0。

表1 25組邊坡的粘聚力和摩擦角取值

?

將25組中的9組不同粘聚力和摩擦角邊坡的破壞特征分析示例結果列于圖4。圖4中No.為對應的邊坡編號，RD分別為計算出的邊坡破壞特征參數，其定義見第二節。從圖4中可以看出，粘聚力和摩擦角的變化對邊坡破壞特征的影響很大。對于編號N=1（即c=5 kPa，φ=10°）的邊坡，邊坡的破壞程度十分明顯，其破壞特征參數的數值都很大。當粘聚力和摩擦角均增加到最大值（即c=25 kPa，φ=30°）時，邊坡的破壞特征參數趨于0，說明具有該土體參數的邊坡處于穩定狀態，其變形很小。

圖4 9組不同粘聚力與摩擦角邊坡的SPH模擬結果

接下來本文分別探究了粘聚力、摩擦角與邊坡的各破壞特征參數RD之間的定量關系。為了便于觀察粘聚力和摩擦角對邊坡破壞參數的影響，本文將粘聚力、摩擦角與對應的邊坡破壞特征參數繪制在同一個三維圖上，如圖5所示。可以看到，隨著粘聚力和摩擦角的增大，邊坡的破壞特征參數RD隨之減小。本文采用線性平面對這些散點進行擬合，得到邊坡粘聚力c、摩擦角φ與邊坡破壞特征參數的關系式，結果如公式（9）所示。

圖5 粘聚力c、摩擦角φ與特征參數RD的關系

研究結果表明：非線性平面可以很好地用于擬合粘聚力、摩擦角和邊坡破壞特征參數的關系，擬合優度R2都高達99%。研究結果表明，對于均質土體邊坡，在本文的計算條件下，單個因素（粘聚力或摩擦角）對邊坡破壞特征參數的影響也是非線性的。該研究結果對于工程實際中的計算和預測邊坡的破壞特征、指導設計、施工和后期加固具均有重要的意義。

4 工程實例

本文選取了廣東省深圳市某排土場作為工程實例進行滑坡模擬。圖6展示了該排土場滑坡前后地形對比圖，該滑坡現場可分為3個區域，由高到低分別為：倒梯形的滑坡物源區、脆弱的滑口區和平緩的堆積區。滑坡物源區原為采石場，工程中期作為排土場傾倒渣土，下方堆積區原為平坦的工業園區，現受到嚴重破壞。缺口出現在該排土場北面，此處為一緩坡，坡高約30 m，滑動土體即余泥渣土從該缺口處涌出，下滑至工業園區，以較高的速度破壞地面建筑物，沖擊原地表土體，并最終形成大面積的堆積體。

圖6 滑坡前后對比圖

以該排土場邊坡作為研究對象，根據勘探資料和該排土場剖面圖，基于SPH法建立動力顯式模型。該排土場模型以平行于邊坡走向為X軸，指向排土場內部為正；垂直邊坡走向為Y軸；豎直方向為Z軸，以向上為正。排土場設計余泥渣土堆積方案為：分10級堆放，每一級高度10 m，坡度1∶2.5，相鄰兩級間設寬3 m的馬道，滑口處于第1級頂面，滑口以下土體未發生滑動，因此可處理為基巖，滑口以上部分全為余泥渣土，高度為90 m，等效坡角為20°。取縱剖面進行分析，其幾何模型如圖7所示。

圖7 邊坡尺寸示意圖（單位：m）

在有限元分析軟件ABAQUS中建立對應的SPH三維邊坡粒子模型，如圖8所示。邊界的大小直接影響應力與應變的分布，本模型的邊界條件為：四周和底部為位移限定邊界，坡面為自由邊界，只允許豎向沉降。網格疏密對計算精度有著明顯影響，網格劃分太稀，則計算誤差大，太密則需要耗費過多的效能，為兼顧數值模擬的效率和計算精度，本文中單元尺寸設為5 m，劃分網格后共49 619個節點，62 318個單元。

圖8 基于SPH粒子的邊坡模型

研究區域的土體自上而下分為余泥渣土、采石場棄土、基巖，坡角處有壓實填土。各土體接觸面狀況復雜，在數值模擬中為了便于建模，將五類不規則土體進行概化，各類土體細觀參數見表2。

表2 不同土體細觀參數取值

圖9為排土場邊坡滑坡后的土體位移云圖，破壞過程共計32 s，滑動土體大部分為原填土臺階上部土體與坡面土體，最大滑動位移為209 m。由完整滑坡模擬過程可知，邊坡的破壞是由坡腳向坡頂逐漸發展，坡面上各點位移沿坡面向上逐漸減小。

圖9 邊坡破壞的SPH粒子位移云圖

5 結論

本文主要針對排土場邊坡土體的粘聚力和摩擦角對均質邊坡破壞特征的變化進行研究，基于拉格朗日型無網格粒子法編寫了數值模擬程序，建立了25組不同粘聚力和摩擦角的SPH邊坡模型和基于SPH法的完整滑坡過程并對其破壞特征參數進行了計算。得出以下結論。

（1）粘聚力和摩擦角的變化對邊坡破壞特征的影響較大。隨著粘聚力和摩擦角的增大，邊坡的破壞特征參數大小隨之降低。

（2）在本文的計算條件下，非線性平面可以很好地用于擬合粘聚力、摩擦角和破壞特征參數RD的關系，擬合優度R2高達99%。

（3）本文的研究結果為定量分析均質邊坡的粘聚力、摩擦角與破壞特征的關系奠定了基礎。為此類排土場邊坡中計算和預測破壞特征、指導設計、施工和后期加固均具有重要的意義。

（4）本文基于排土場的邊坡實況，建立邊坡的非線性有限元模型，數值計算結果直觀地表現了邊坡的變形破壞情況和土體內部應力、應變狀態，再現了滑坡從穩態到失穩破壞的全過程。

（5）本文的研究對象是三維均質的理想模型，然而自然界中邊坡是非均質的，且邊坡形狀、土體分布各異，這也是本文研究中的一個不足之處。在后續研究中，將繼續深入研究真實形態特征的邊坡，并考慮其土體力學參數的空間變異性與其邊坡破壞特征的關系。