基于顆粒流方法的滑坡災變過程及沖擊致災預測研究

陳 鑫，趙 洲，楊 鵬

(1.西安科技大學地質與環境學院，陜西 西安 710054；2.西安地質礦產勘查開發院有限公司，陜西 西安 710100；3.中國有色金屬工業西安勘察設計研究院有限公司，陜西 西安 710054)

0 引 言

滑坡是指由于受到地層巖性、水、地震及人類工程活動等因素影響，坡體沿貫通的剪切破壞面或帶，以一定加速度下滑移動的不良地質現象[1]。隨著我國經濟的發展及各項基礎設施建設向城市臨近的山區不斷拓展，滑坡問題越顯突出，對人類生命和財產安全造成極大危害。據統計，全國平均每年地質災害直接經濟損失達200多億元，其中滑坡造成的經濟損失達2/3。因此，滑坡研究和防治工作一直以來是國家關注的重點。近年來，關于滑坡破壞運動過程及沖擊致災研究受到許多國內外專家學者的普遍關注[2- 6]。黃潤秋等[7]采用有限元數值模擬方法，對膨脹土滑坡的特點和發生機理分別從強度分層、裂隙影響和蠕變效應進行的研究得出，膨脹土路塹邊坡滑坡受上述3種因素共同作用，具體表現為淺層、小規模、漸進性和滑坡后緣通常出現在坡腰部位。亓星等[8]通過對甘肅黑方臺滑坡野外調查和位移監測，探討了滑坡變形破壞過程，并對滑坡的發生機制進行了詳細解釋。吳越等[9]利用光傳感器和計時器組合而成的摩擦系數測量系統，根據相似原理建立滑坡滑體下滑沖擊模型，探討了滑體下滑和沖擊致災全過程中的能量耗散規律，為準確計算滑體沖擊能奠定了基礎。薛強等[10]采用三維數字高程，提出了滑坡強度定量評估技術方法，并對黑方臺南緣區32處滑坡強度進行了評估，為滑坡風險評估提供了新思路。黨杰等[11]根據物理學中內能轉化原理，以地質災害發生時釋放的能量作為強度評估的依據，建立了滑坡強度與能量關系，提出了單體滑坡強度評估的指標和方法。以上滑坡研究方法大多處于定性或半定量研究，而定量研究方法仍然處于緩慢發展階段。

“5·12”地震后，受地震影響，樓子溝滑坡土體產生松動變形，地面隆起，房屋裂縫變形加劇，少數房屋即將倒塌，滑坡后部出現弧形裂縫，醉漢林密布，嚴重影響21戶、90人、100余間房及百余畝良田的安全，可造成經濟損失達400萬元以上。為此，本文采用PFC2D數值模擬方法，以樓子溝滑坡為研究對象，仿真動態模擬滑坡發生全過程，對滑坡沖擊速度、沖擊距離、沖擊力大小等指標進行定量分析，同時對滑坡沖擊致災程度進行定量化評價。

1 滑坡概況

樓子溝滑坡位于陜西省漢中市勉縣老道寺鎮樓子溝村一組，屬北部山區南緣丘陵區。滑坡平面呈一定弧長的長條形，主滑方向90°，滑坡體長約120 m，前緣寬度約360 m，后緣寬度約300 m，平均厚度約5 m，體積約2.16×105m3，屬中型膨脹土滑坡。滑坡后緣以陡坡為界，陡坎高度約5～8 m；滑坡左、右兩側以小沖溝為界；滑坡前緣為村民水田，在邊界外側未發現變形跡象，為原始地貌；滑坡中后部地表可見明顯下挫，邊界清晰；滑坡后緣有泉水出露，且有一灌溉渠，渠中常年有水，水量較小。滑坡平面見圖1。

根據收集已有資料和勘查結果，地層自上而下為第四紀堆積粉質粘土、下更新統沖洪積粉質粘土、第三系泥巖。滑體物質主要為滑坡堆積粉質粘土，黃褐色，土體結構松散，節理裂隙十分發育，土質呈斑塊狀，掰開后形如石榴粒，可見明顯的灰褐色膨脹土礦物，網格狀裂縫十分發育，最大寬度3 mm，灰色粘土薄膜附于裂隙表面，表面擦痕明顯，裂隙面光滑，手捻有滑感。滑帶土由灰褐色粉質粘土組成，呈可塑狀，滑膩感很強，滑面光滑，可見擦痕。

圖1 滑坡平面

2 滑坡模型構建

2.1 滑坡巖土體宏細觀力學參數標定

顆粒流(PFC)方法采用中心差分法，把虎克定律和牛頓第二定律遍歷整個顆粒集合體，模擬顆粒間的運動及相互作用，運用軟件特定的粘結接觸模型將其集合起來模擬巖土體的完整性，不受變形量的限制，可以方便地處理非連續介質力學問題，有效地模擬介質開裂、分離等現象[12]。在實際工程中，通過室內物理試驗得到的宏觀力學參數(如變形模量、內摩擦角、粘聚力等)可以直接在工程實際中運用計算，而在顆粒流數值模擬方法中，巖土體材料的宏觀力學性質是通過顆粒總體的運動軌跡來反映的，而其運動受顆粒細觀力學參數(如顆粒接觸剛度、顆粒粘結強度、顆粒粘結剛度等)的影響。因此，細觀力學參數的確定是顆粒流數值模擬分析的關鍵。鑒于關于顆粒細觀力學參數和材料宏觀力學參數之間的定量關系還沒有確定的公式[13]，故PDC2D中需要采用雙軸壓縮試驗模擬室內物理試驗，以此得到細觀力學參數與宏觀力學參數之間的對應關系[14]。

本文采用的雙軸壓縮試驗模型見圖2。模型尺寸為10 m×5 m，顆粒粒徑為0.1～0.15 m，生成顆粒總數為885個。滑坡巖土體顆粒細觀力學參數見表1。在表1的細觀力學參數下，施加圍壓分別為50、100 kPa和150 kPa，得到應力-應變關系曲線(見圖3)，再根據M-C理論繪制M-C包絡線(見圖4)，最后通過計算得到雙軸壓縮試驗中材料的宏觀力學參數。綜上，通過數值模擬得到的巖土體宏觀力學參數粘聚力c=13.3 kPa，內摩擦角φ=12°，與室內物理試驗中c=11.5～14.7 kPa，φ=11.0°～15.3°對比可知，數值模擬所得細觀力學參數符合要求。

表1 滑坡巖土體顆粒細觀力學參數

圖2 雙軸試驗模型(單位：m)

圖3 應力-應變關系

圖4 M-C包絡線

2.2 滑坡PFC2D模型構建

PFC中提供了10種內嵌的接觸模型，本次數值模擬采用平行粘結模型，因為其既可以傳遞顆粒間的作用力，也可以傳遞顆粒間的作用力矩，從而更好的模擬顆粒間粘結破壞。本文采用PFC中ball-ball方法建立滑坡模型[15]，具體步驟如下：①根據滑坡主剖面圖，按1∶1在PFC中建立矩形區域并生成顆粒；②通過命令導入AutoCAD中建立的滑坡剖面圖并生成滑坡邊界墻體使顆粒分組(見圖5)；③利用PFC中ball delete range geometry命令刪除削坡部位顆粒；④根據表1中顆粒細觀力學參數對顆粒進行賦值，同時對顆粒施加重力作用使滑坡在重力作用下達到初始應力平衡狀態。

圖5 未削坡前滑坡模型

根據上述步驟得到樓子溝滑坡PFC2D模型，見圖6。模型長165 m，高18.5 m，顆粒總數28 270個。其中，滑體顆粒總數7 058個，滑床顆粒總數21 212個。

圖6 滑坡PFC2D模型

3 滑坡災變過程及沖擊致災預測

3.1 滑坡災變過程分析

為了更好地分析樓子溝滑坡變形特征及運動過程，根據滑坡模型尺寸及測量圓設置原則[16]，共建立4個測量圓，以此來監測滑坡滑體部位應力變化情況。滑坡位移見圖7。從圖7可知，當運行至4 500步時，滑坡坡腳出現剪切破壞，滑坡后緣出現拉張破壞，這主要是由于降雨導致滑體自重不斷增加，當其超過鎖固段強度時，滑坡發生剪出破壞；當運行至180 000步時，由于顆粒間相互摩擦碰撞消耗能量，滑坡動能逐漸減少，最終停止運動堆積在坡腳處。通過測量，模擬中顆粒最大位移為58.8 m，滑坡實際位移為36.9 m。

圖7 滑坡位移

滑坡開始發生時，坡腳處剪切破壞。應力變化見圖8。從圖8可知， 0～13 s內1號測量圓中應力持續增大，由于滑坡體坡度較緩，故x方向的應力增量大于y方向，且近似為2倍；20 s后滑坡運動停止，應力達到平衡。12 s時2號測量圓內x方向應力達到峰值110 kPa，約為y方向應力的2.5倍，此時y方向應力來源主要為上部滑體下滑堆積擠壓產生；3號和4號測量圓內應力逐漸減小，最終在11 s和4 s時分別歸0，這是由于滑坡整體發生破壞，3號和4號測量圓位置處滑體部位由滑床滑出，故應力為0。

圖8 應力變化

滑坡速度及位移分析是研究滑坡破壞運動過程的重要途徑。樓子溝滑坡全長約120 m，本文以30 m為1段監測部位，每個部位分上、中、下3排，每排設置監測顆粒5個，對滑坡速度及位移進行定量描述，以加深對樓子溝滑坡破壞運動過程的認識。滑坡監測結果見圖9。從圖9可知：

圖9 滑坡監測結果

(1)各個監測點運動主要分為加速階段、過渡階段、減速階段。滑坡開始破壞時發生整體同等加速運動，滑坡體0～30 m內顆粒最大速度為3.6 m/s，隨著滑坡前端剪出破壞，30～120 m范圍內滑體由于前部失去支撐而向下運動，顆粒速度逐漸增大，最大速度達6.6 m/s。滑坡體0～30 m監測點顆粒減速階段的歷時遠大于其他范圍內監測點，這是由于前部顆粒不僅受30～90 m范圍內顆粒的推動獲得能量，同時還傳遞90～120 m范圍內顆粒的能量，而90～120 m范圍內顆粒由于沒有后續能量補給，無法維持遠距離運動，故在減速階段迅速停止運動。

(2)滑坡體0～30 m范圍內監測點顆粒最大位移為34.8 m；30～120 m范圍內顆粒位移呈現依次增大特點，最大顆粒位移為58.8 m。該運動特征符合滑坡下部剪出破壞，逐級牽引滑坡體上部向下滑動的復合式滑坡規律。

綜上所述，樓子溝滑坡運動過程可分為0～4 s等速變形階段、4～15 s加速變形階段、15 s之后堆積靜止階段3個階段。滑坡破壞初期主要以蠕滑變形為主，隨著形變量的增加，滑坡鎖固段破壞，滑坡下部剪出，逐級牽引上部滑體順破裂面下滑，最終堆積于坡腳處，滑坡全過程符合典型牽引式滑坡的運動特征。

3.2 滑坡沖擊致災預測

滑坡沖擊致災研究是滑坡風險評價的關鍵環節[17]，沖擊力的大小直接反映了建筑物等承載體損毀的嚴重程度。因此，本文通過在滑坡不同坡腳距位置設置剛性擋墻代替建筑物，建立滑坡沖擊模型(見圖10)，以此統計滑坡下滑時對擋墻的沖擊力，為滑坡沖擊致災分析提供依據。

圖10 滑坡沖擊模型

滑坡沖擊力見圖11。從圖11可知，滑坡沖擊擋墻過程可分為3個階段：①自由下滑階段。滑坡啟動但未接觸擋墻。②沖擊階段。滑坡不斷沖擊擋墻，滑體顆粒間及顆粒于擋墻相互碰撞耗能使滑坡能量逐漸減小。③沖擊靜止階段。此時滑坡運動基本停止，滑體堆積于坡腳和擋墻之間，沖擊力來自滑體堆積后對擋墻的靜止土壓力，為恒定值。

圖11 滑坡典型沖擊力

滑坡在沖擊階段內的沖量E可以根據沖擊曲線和時間軸的面積得出，再根據等效力在相同時間內沖量相等的原理求出等效沖擊力F，具體關系如下

(1)

對滑坡破壞全過程模擬分析可知，滑坡最大位移為36.9 m，37 m處沖擊力為0。故從坡腳開始，每隔2 m設置剛性擋墻，共建立滑坡沖擊模型19個，以此監測統計滑坡下滑對不同坡腳距處擋墻的沖擊力。利用Origin數據分析軟件計算得到不同坡腳距處等效沖擊力數據，最后線性擬合繪制出滑坡坡腳距-沖擊力關系(見圖12)。由圖12可知，滑坡等效沖擊力隨坡腳距的增大而逐漸遞減，具體關系為F=-15 188.608 96x+543 387.274 59(R2=0.969 99)。滑坡對建筑物等承災體的損毀程度隨坡腳距增大而逐漸減弱。

圖12 滑坡坡腳距-沖擊力關系

根據現場實際調查資料，滑坡區內房屋大多以磚混結構為主，由于砌體結構主要發生脆性變形，故建筑物等承災體不滿足設計要求時認為是失效狀態，等效沖擊力作用位置由滑坡沖擊結束時沖擊高度及墻體所受荷載分布形式共同決定，由此建立建筑物失效模型，參考吳越等[18]承災體易損性研究可定義滑坡沖擊致災指標R，具體關系如下

(2)

式中，F′為建筑物所能承受的極限荷載值。

由GB 50003—2011《砌體結構設計規范》可知，受彎構件承載力及受剪承載力，應按下列公式計算

M≤ftmW

(3)

V≤fvbz

(4)

Z=I/S

(5)

式中，M為彎矩設計值；ftm為砌體彎曲抗拉強度設計值；W為截面抵抗矩；V為剪力設計值；fv為砌體結構抗剪強度設計值；b為截面寬度；Z為內力臂；I為截面慣性矩；S為截面面積矩。

本文取單面墻寬，分別計算墻體受彎矩控制的荷載值F彎和剪力控制的荷載值F剪并取其中最小值作為建筑物所能承受的極限荷載值F′，即F′=min(F彎，F剪)。設置墻體距離坡腳20 m，根據圖11方法計算得到此時滑坡等效沖擊力F=2.45×105N，根據勘察資料，墻體材料為MU10普通燒結磚和M7.5水泥混合砂漿，墻體采用“三七墻”方式砌筑，墻體厚度365 mm，由式(3)、(4)、(5)計算得到墻體所能承受極限荷載值F′=min(F彎，F剪)=F彎=2.04×103N，將F、F′代入式(2)計算得到滑坡沖擊致災指標R=120.1。根據致災程度分析可知，滑坡發生時墻體無法承受沖擊損毀，會造成嚴重后果，對該滑坡應進行加固治理。

4 結 語

本文以漢中市勉縣老道寺鎮樓子溝滑坡為研究對象，利用PFC2D數值模擬方法對滑坡破壞運動災變全過程進行了仿真模擬，得出以下結論：

(1)滑坡破壞初期主要以蠕滑變形為主，隨著形變量的增加，滑坡鎖固段破壞，逐級牽引上部滑體順破裂面下滑，表現為典型牽引式滑坡的運動特征。

(2)滑坡體最大致災位移為36.9 m，最大運動速度為6.6 m/s，等效沖擊力隨坡腳距增大呈線性遞減特征，該特征可為滑坡沖擊致災定量研究提供參考。

(3)在滑坡致災范圍內，坡腳距離為20 m時，滑坡沖擊致災指標R=120.1，滑坡變形破壞對下方建筑物必然造成損毀，應對滑坡體進行加固治理，防止滑坡下滑造成嚴重的人員傷亡事故。