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

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

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

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

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

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

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

1 SPH法

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

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

在SPH中，計算域使用粒子進行劃分，粒子攜帶場變量并隨材料移動?？刂品匠炭梢酝ㄟ^跟蹤粒子的運動和攜帶變量的變化來求解。要做到這一點，必須使用SPH插值離散化控制方程。場函數(shù)可以由以下積分插值近似：

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

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

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

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

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

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

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

2 模型構(gòu)建與邊坡破壞特征量化方法

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

圖1 邊坡尺寸示意圖

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

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

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

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

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

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

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

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

?

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

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

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

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

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

4 工程實例

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

圖6 滑坡前后對比圖

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

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

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

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

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

表2 不同土體細觀參數(shù)取值

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

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

5 結(jié)論

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

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

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

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

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

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