魏良針，王建平，楊 磊，應(yīng)濤濤，龔禮岳*，吳則祥

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

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

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

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

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

1 SPH法

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

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

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

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

類似的，場(chǎng)函數(shù)的空間導(dǎo)數(shù)可以通過(guò)以下方程近似：

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

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

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

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

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

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

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

圖1 邊坡尺寸示意圖

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

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

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

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

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

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

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

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

?

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

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

接下來(lái)本文分別探究了粘聚力、摩擦角與邊坡的各破壞特征參數(shù)RD之間的定量關(guān)系。為了便于觀察粘聚力和摩擦角對(duì)邊坡破壞參數(shù)的影響，本文將粘聚力、摩擦角與對(duì)應(yīng)的邊坡破壞特征參數(shù)繪制在同一個(gè)三維圖上，如圖5所示?？梢钥吹剑S著粘聚力和摩擦角的增大，邊坡的破壞特征參數(shù)RD隨之減小。本文采用線性平面對(duì)這些散點(diǎn)進(jìn)行擬合，得到邊坡粘聚力c、摩擦角φ與邊坡破壞特征參數(shù)的關(guān)系式，結(jié)果如公式（9）所示。

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

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

4 工程實(shí)例

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

圖6 滑坡前后對(duì)比圖

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

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