朱志男，李慧亭，陳 羽，劉太輝，成昊凌

(1.中國建筑第八工程局有限公司，甘肅 蘭州 730030；2.西安交通大學 人居環(huán)境與建筑工程學院，陜西 西安 710049)

滑坡災害對人類的生命財產(chǎn)安全造成了極大的威脅，根據(jù)國內外的滑坡記錄，許多滑坡都是由降雨導致[1]。2011年“9.16”特大暴雨誘發(fā)了四川省南江縣數(shù)以千計的群發(fā)性緩傾角淺層土質滑坡[2]，2013年7月，在持續(xù)二十多天的降雨作用下，陜西省延安市全市有7 594處發(fā)生滑塌，造成了42人遇難，經(jīng)濟損失達66.15億元[3-4]。可見，研究降雨誘發(fā)滑坡的機理、分析邊坡滑坡過程對提出邊坡防災減災的對策具有重要意義。

目前國內外學者對降雨誘發(fā)滑坡做了大量的研究，研究方法主要為理論分析、室內試驗與數(shù)值模擬[5]，由于數(shù)值模擬成本低、效果顯著的特點被廣泛應用[6]。許方領等通過有限元軟件ABAQUS軟件分析了降雨入滲后的滲流場及位移和安全系數(shù)的變化規(guī)律[7]，王磊等利用有限元軟件Geostudio軟件分析了降雨條件下坡體含水率和孔壓響應情況，并探討了裂縫對陡坡穩(wěn)定性的影響[8]。通過有限元方法模擬滲流比較簡單高效，但有限元方法難以考慮大變形，無法反映滑坡過程中裂隙的延伸過程，可以通過顆粒流方法來分析邊坡的滑坡過程，近年來眾多學者開始利用顆粒離散元的方法來研究邊坡的穩(wěn)定性，汪華安通過強度折減的方法借助顆粒流數(shù)值模擬研究了煙墩嶺邊坡在降雨工況下的變形破壞模式[9]，張家勇以貴陽省開陽縣魚鰍坡滑坡為研究對象，通過PFC3D分析了滑坡的破壞過程[10]；吳謙通過PFC3D對黃土邊坡降雨沖刷過程進行了流固耦合模擬，研究了降雨過程中邊坡遭受侵蝕程度及水流侵蝕能力的分布規(guī)律[11]；但顆粒流方法缺乏有效的本構模型來反應飽和-非飽和應力應變關系，現(xiàn)有的PFC-CFD耦合方法計算效率較低，且流體網(wǎng)格形狀較為規(guī)則單調，難以應用于工程實際。

為解決有限元方法無法計算大變形以及顆粒流方法中流體-顆粒耦合計算較為復雜的問題，本文將有限元與離散元相結合，提出了一種PFC-Geostudio耦合的方法，將有限元軟件Geostudio的滲流計算結果導入到顆粒離散元方法中，通過等效方法來研究邊坡在降雨條件下的變形情況以及滑坡過程。

1 PFC-Geostudio邊坡降雨模擬方法

降雨對邊坡土體有三個方面的影響：(1)土體由天然狀態(tài)轉化為飽和/非飽和狀態(tài)，重度增加，增大了下滑力；(2)土體強度參數(shù)減小，如摩擦角、內聚力等，增大滑坡風險；(3)土體內發(fā)生滲流，土顆粒受到滲流力的作用，當降雨強度較大時，坡表還會發(fā)生徑流，對土顆粒產(chǎn)生拖曳力作用[12]。因此，可通過Geostudio來計算滲流場，然后將滲流計算的結果導入至PFC方法中，根據(jù)顆粒所在位置的含水率，改變其重度與接觸參數(shù)，實現(xiàn)降雨導致土體重度和強度參數(shù)的變化模擬；對滲流區(qū)域施加滲流力，當降雨強度較大時，額外考慮拖曳力，以實現(xiàn)降雨對土體的力學作用的模擬；通過以上方法即可在顆粒流中等效實現(xiàn)降雨效果，隨后可進行邊坡變形與滑坡分析，具體流程可見圖1。

圖1 Geostudio滲流計算模型圖Fig.1 Geostudio-PFC calculation method

1.1 Geostudio滲流計算

有限元軟件Geostudio內置的SEEP/W模塊，通過非飽和土力學計算方法，可以快速模擬邊坡在降雨條件下滲流場變化情況，并得到邊坡內浸潤線分布以及各點的含水量。為研究在降雨條件下邊坡內的滲流場，建立如圖 2所示的有限元模型，模型長度為50 m，高度為20 m，單元網(wǎng)格長度設置為0.3 m，整個邊坡模型劃分為8 676個單元，8 879個節(jié)點，將模型的上邊界設置為降雨邊界，兩側及底部設置為零流量邊界，滲流計算中所使用的參數(shù)如表1所示。

圖2 Geostudio滲流計算模型圖Fig.2 Geostudio seepage calculation model diagram

根據(jù)氣象部門對降雨強度的劃分，降雨通?？煞譃闉樾∮?、中雨、大雨、暴雨、大暴雨和特大暴雨六個等級，為研究在降雨條件下邊坡的滲流場情況，設定降雨強度為30 mm/24 h(大雨)，原始水位為地下30 m，連續(xù)48 h降雨后，邊坡滲流場如圖3(a)所示，可見沿雨水入滲方向，形成了一定厚度的飽和帶，繼續(xù)向下為非飽和帶，最下側土體為天然狀態(tài)，未受雨水入滲影響，提取邊坡前緣、邊坡底部、邊坡中部、邊坡頂部、邊坡后緣五個位置沿入滲方向的各點的含水率，形成圖3(b)所示曲線圖，可見對于圖2所示的均質土體邊坡，降雨入滲規(guī)律與土體的所在位置無關，僅與其所在深度相關，土體含水率最大的區(qū)域出現(xiàn)在飽和帶，最大含水率為40%，飽和帶的寬度約為1.6 m，在非飽和帶，土體含水率沿入滲方向持續(xù)減小，非飽和帶的厚度約為1.8 m，含水率最小為15.5%。

表1 Geostudio中的模型參數(shù)Tab.1 Model parameters in Geostudio

圖3 降雨影響下邊坡滲流場情況Fig.3 Slope seepage field under the influence of rainfall

1.2 PFC-Geostudio耦合方法

為實現(xiàn)PFC-Geostudio的耦合，可在PFC中建立相同尺度的模型，將Geostudio中的計算結果，即每個位置的含水率情況導入；在顆粒流模型中，根據(jù)各個顆粒的位置，修改其重度，若顆粒處于飽和區(qū)域，則顆粒的重度應取飽和重度，顆粒的接觸參數(shù)也需按照飽和參數(shù)進行折減，根據(jù)文獻[13]的統(tǒng)計結果，飽和參數(shù)可按照天然參數(shù)的15%折減，若顆粒處于非飽和區(qū)域，顆粒的重度應根據(jù)其所在位置的含水率進行折減，其接觸參數(shù)也應進行相應的折減，具體可見下式：

(1)

(2)

式中：γ1—飽和區(qū)內顆粒重度，N/m3；γ2—非飽和區(qū)顆粒重度，N/m3；γ3—天然區(qū)域顆粒重度，N/m3；ω—顆粒所在位置的含水率，可從Geostudio中的計算結果獲?。沪?—飽和含水率；ω3—天然含水率；f1—飽和區(qū)顆粒的接觸參數(shù)(此處以摩擦系數(shù)為例)；f3—天然區(qū)域顆粒的接觸參數(shù)。

除改變顆粒重度與接觸參數(shù)外，在飽和區(qū)域由于雨水入滲，顆粒還會受到滲流力的作用，顆粒所受到的滲流力可簡化等效為[14-16]:

J=γfiV

(3)

式中：J—滲流力，N；方向為沿雨水入滲方向；γf—流體重度，N/m3；i—流體的水力梯度，V—流體體積，m3。

當降雨強度較大，超過土體的入滲能力，則在邊坡表面會發(fā)生徑流。在水的沖刷作用下，坡表土體會受到?jīng)_刷作用，則土顆粒會受到拖曳力的作用，其方向為沿雨水沖刷方向，其大小為：

fdrag=f0ε-χ

(4)

式中：f0—單個顆粒所受的拖曳力，N；ε—顆粒所在流體單元的孔隙率，ε-χ為考慮局部孔隙度的經(jīng)驗系數(shù)。這個修正項可以使拖曳力同時適用于孔隙度較高和孔隙度較低的系統(tǒng)，而且流體的雷諾數(shù)也可以大幅度取值。單個圓形顆粒所受流體拖曳力為[17]：

(5)

式中：Cd—拖曳力系數(shù)；ρf—流體密度，kg/ m3；r為顆粒半徑，m；u—顆粒速度矢量，m/s；v—流體速度矢量，m/s；不同位置的坡面徑流速度各不相同，當水沿坡面流下時，由于重力的作用，水流的速度從坡頂?shù)狡碌壮掷m(xù)增加，在邊坡底部水流速度達到最大，當水流到達坡底時，由于土顆粒的阻礙，水在坡底流動時，其速度又會持續(xù)減小。

如圖4所示，位于坡面的顆粒A處的徑流速度大小與方向為:

圖4 不同位置坡面徑流速度示意圖Fig.4 Schematic diagram of the velocity of slope runoff at different positions

(6)

(7)

(8)

式中：vA—顆粒A處的徑流速度，m/s；vmax—坡表徑流的最大沖刷速度，m/s；如圖 4所示，h1—顆粒A距坡底的高度，m；h2—坡頂距坡底的高度，m；vAx—顆粒A處徑流速度沿x方向的分量，m/s；vAy—顆粒A處徑流速度沿y方向的分量，m/s；Δx—水平方向的單位增量;Δy—水平方向增量為Δx時坡面線的高程增量。

位于坡底的顆粒B處的徑流速度大小為:

(9)

式中：vB—顆粒B處的徑流速度，方向為沿徑流的沖刷方向，m/s；x1—顆粒B距坡底的水平距離，m；x2—坡底的長度，m。

式(5)中，拖曳力系數(shù)Cd為：

(10)

式(5)中，經(jīng)驗系數(shù)χ采用下式進行計算：

(11)

式(10)與式(11)中，Rep表示顆粒雷諾數(shù)，并采用下式進行計算：

(12)

式中，μ—流體動力粘度，N·s/m2。

2 模型建立與實現(xiàn)方法

2.1 顆粒流模型建立及參數(shù)標定

在顆粒流方法中，首先導入邊坡的邊界生成墻體，然后在邊界內生成完整的細觀顆粒體系。當按照真實土顆粒的粒徑生成模型時，會導致顆粒數(shù)量過多，嚴重影響計算效率，但如果顆粒粒徑取值過大，生成顆粒太少，模型的精度也會受到影響，因此建立顆粒流模型時，選擇一個合適的粒徑范圍至關重要，根據(jù)孫其誠[18]的研究，當模型的特征長度比值L/R≥50時，顆粒粒徑對材料的力學性質影響較小，因此選取顆粒粒徑為0.04～0.08 m，共生成顆粒83 771個，由于通過初始命令生成的模型比很不均勻，會出現(xiàn)部分區(qū)域顆粒之間空隙較大、部分區(qū)域顆粒之間應力集中的現(xiàn)象[19]，因此需要運用伺服機制讓模型變的更加均勻，伺服后的模型及其力鏈分布情況如圖5所示。

圖5 顆粒流初始模型Fig.5 Initial model of particle flow

在顆粒流方法中，通過顆粒(ball)來模擬材料顆粒，通過接觸(contact)來模擬顆粒之間的相互作用，通過賦予顆粒不同的細觀參數(shù)來模擬顆粒之間不同的力學性質，因此選擇合適的接觸模型以及對接觸賦予合適的參數(shù)至關重要[20]。

當采用線性接觸模型，顆粒間的接觸處于粘結狀態(tài)時，接觸模型是線彈性的，直到受力狀態(tài)超出強度極限以及粘結被破壞才會使界面處于非粘結狀態(tài)，在非粘結狀態(tài)下，線性接觸模型處于線彈性和摩擦狀態(tài)，當剪切力超出庫倫極限時就會發(fā)生滑動，該粘結模型能很好的模擬土體應力應變特性及其破壞特性，因此顆粒間的接觸采用線性接觸模型。

由于賦予顆粒的細觀參數(shù)與材料的宏觀參數(shù)之間的關系并不是一一對應的，因此為確保顆粒流模型細觀參數(shù)的合理性，需要進行參數(shù)標定，即不斷調試模型的細觀參數(shù)，直到模型的宏觀響應與我們所需的模型宏觀參數(shù)相接近時，參數(shù)標定方才結束。構建如圖 6(a)所示的巖土試樣，試樣的孔隙率和粒徑取值與邊坡模型取值相同，粒徑取0.04～0.08 m，孔隙率取0.17。對建立的巖土試樣模型進行伺服，使試樣模型應力分布均勻，隨后賦予接觸參數(shù)，對巖土試樣進行不同圍壓下的雙軸壓縮試驗，提取模型的應力應變曲線(圖 6(a))，繪制莫爾應力圓，如圖 6(b)所示，即可得試樣的宏觀參數(shù)，不斷調整接觸的細觀參數(shù)，使得試樣的宏觀參數(shù)為c=23 kPa，φ=14 °，邊坡模型的細觀參數(shù)最終取值如表2所示。

圖6 參數(shù)標定圖示Fig.6 Parameter calibration diagram

表2 顆粒流模型細觀力學參數(shù)Tab.2 The meso-mechanical parameters of the calculation model

2.2 顆粒流模擬中降雨實現(xiàn)方法

經(jīng)參數(shù)標定后，對建立的初始邊坡模型賦予參數(shù)，施加重力，運行至平衡，隨后對所有顆粒的速度和位移清零，此時的模型即為天然狀態(tài)下的邊坡模型。設定降雨強度為30 mm/24 h，經(jīng)Geostudio滲流計算，將各點的含水率導入顆粒流模型中，遍歷所有顆粒，根據(jù)顆粒所在位置的含水率按照式(1)和式(2)修改顆粒的重度與接觸參數(shù)，實現(xiàn)效果如圖7所示。

當對飽和區(qū)域的顆粒施加滲流力時，首先應篩選出受滲流力影響的顆粒。遍歷所有顆粒，由于飽和土的含水率為40%，當顆粒所在位置的含水率大于39%時，即可判定其屬于飽和區(qū)，會受到滲流力的作用，即可按照式(3)對顆粒施加等效滲流力。當降雨強度大于土體的入滲能力時，坡表由于雨水沖刷，坡表顆粒會受到拖曳力的作用，因此在具體施加拖曳力時，首先需要判斷降雨強度是否大于土體入滲能力，僅當降雨強度大于土體入滲能力時才調用拖曳力施加函數(shù)，在進行拖曳力的具體施加時，遍歷所有顆粒，判定出坡表受徑流影響的顆粒，如圖8(a)所示，由于拖曳力的方向為雨水徑流方向，因此當顆粒位于坡底時，徑流的方向為沿水平方向，顆粒位于坡面時，徑流的方向為沿坡面方向，不同位置的徑流速度的大小與方向可按照圖4以及式(6)—(9)考慮，再通過Fish函數(shù)獲取顆粒的粒徑、速度、位置等信息，隨后即可根據(jù)式(4)—(12)來施加拖曳力。當降雨強度為100 mm/24 h時，顆粒流方法中滲流力與拖曳力的施加效果如圖8(b)所示。

3 結果分析

3.1 大雨工況下邊坡變形分析

為研究在降雨條件下邊坡的變形情況，設定降雨強度為30 mm/24 h，降雨持續(xù)時間為48 h，經(jīng)Geostudio滲流計算，計算結果如圖3(a)所示，將各個位置的含水率導入PFC中，通過1.2節(jié)所述的耦合方法來等效實現(xiàn)降雨效果。顆粒流邊坡在降雨條件下的變形情況如圖9(a)所示，可見最大位移區(qū)域位于坡頂，最大位移為8.3 cm，變形主要集中在坡頂飽和區(qū)域，土體位移方向如圖9(b)所示，可見位移方向主要為土體沉降，產(chǎn)生沉降的原因在于降雨導致土體容重增加以及土體強度參數(shù)的減弱，此外坡面土體由于左側無土體支撐，因此會產(chǎn)生一定沿坡面滑動的趨勢，若降雨強度持續(xù)增加，則降雨入滲深度會繼續(xù)增大，同時飽和帶與非飽和帶的厚度也會持續(xù)增加，會導致坡面土體穩(wěn)定性會進一步降低，有潛在滑坡的趨勢。

圖10 降雨條件下的測點位移Fig.10 Measuring points displacement under rainfall conditions

圖11 暴雨工況下邊坡滑坡過程及裂隙延伸情況Fig.11 Landslide process and fissure extension of the slope under rainstorm conditions

在坡底、坡面、以及坡頂設置監(jiān)測點，以進一步了解在計算過程中邊坡不同位置土體的變形情況。由于單獨監(jiān)測一個顆粒隨機性較強，可能會出現(xiàn)少數(shù)顆粒飛出模型邊界的情況，因此選擇數(shù)個顆粒作為一個測點，形成一個顆粒團，取其平均位移以反映測點所在位置土體的位移，可更準確的反映土體的變形情況，且當邊坡發(fā)生失穩(wěn)破壞時，測點顆?？呻S之運動，依舊能反映測點所在區(qū)域邊坡的變形情況，如圖10(a)所示，沿坡面共設置測點15個，每個測點半徑為0.3 m，平均每個測點包含顆粒20個。當降雨強度為30 mm/24 h、降雨持續(xù)時間為48 h時，經(jīng)PFC-Geostudio耦合計算后，提取各測點位移，形成如圖10 (b)所示曲線圖，可見沿水平方向，坡底和坡頂各測點位移變化不大，但坡頂位移整體大于坡底，坡面各測點的位移變化較大，各測點的位移隨高程增加而增加，最大位移出現(xiàn)在坡頂位置。

3.2 暴雨工況下邊坡滑坡過程分析

為研究邊坡在暴雨工況下的變形破壞情況，設定降雨強度為100 mm/24 h，降雨持續(xù)時間為48 h，此時由于降雨強度遠大于土體的入滲速度，因此在坡面會產(chǎn)生徑流，需額外考慮拖曳力對坡表顆粒的影響，在暴雨工況下，拖曳力的施加效果如圖8(b)所示，經(jīng)PFC-Geostudio耦合計算，邊坡的變形破壞過程如圖11所示，邊坡的滑動過程可劃分為3個階段。

第一階段如圖11(a)—(b)所示，坡面顆粒由于雨水沖刷作用，在坡腳首先發(fā)生破壞，底部顆粒在拖曳力的作用下脫離邊坡，沿坡面滑動至坡底，同時邊坡內飽和區(qū)域形成大量裂隙，但并未滑動；

第二階段如圖11(b)—(d)所示，由于坡腳顆粒沖刷滑離邊坡，導致上部土體失去支撐，首先在坡腳發(fā)生小范圍滑坡，隨后裂隙拓展延伸，從坡腳已滑動區(qū)域的頂部延伸至坡頂，形成大范圍滑坡，在滑動過程中，滑體內部出現(xiàn)大量裂隙，滑體逐漸破碎為散體；

第三階段如圖11(d)—(g)所示，邊坡后緣土體由于失去前側土體的支撐，內部出現(xiàn)大量裂隙，在自重作用下失穩(wěn)，出現(xiàn)了多次小規(guī)模的滑動，邊坡滑坡的最終狀態(tài)如圖11(h)所示，坡面土體大范圍滑動，滑體破碎并堆積至坡底。

為確定在暴雨條件下坡面不同位置土體的滑坡時序，提取在滑坡起始階段中測點6、測點8、測點10的速度信息，如圖12所示，可見由于坡面徑流沖刷作用的影響，位于坡腳的測點6速度最先增加，且速度較大，當坡腳顆?；瑒雍?，上部土體逐步變形，位于坡面中部的測點8的速度開始逐步增加，隨后測點10的速度隨之增加，但其速度小于測點6，可見測點8附近土體先于測點10附近土體滑動，滑坡過程為明顯的漸進式牽引式滑坡。

圖12 滑坡過程各測點位移Fig.12 Displacement of each measuring point during the landslide process

文獻[21]搭建室內降雨滑坡框架模型，通過在模型頂部設置噴頭來模擬降雨，試驗中的滑坡過程如圖13(a)—(d)所示，在降雨沖刷作用下表面土體流失，從坡腳處開始發(fā)生破壞，產(chǎn)生滑坡臨空面，逐步誘發(fā)了大規(guī)模滑坡，并引發(fā)了數(shù)次小規(guī)模滑坡，滑坡過程為典型的牽引式滑坡，與上述采用耦合計算方法得到的滑坡過程(圖13(e)—(h))基本一致，可證實本方法的可行性。

圖13 滑坡過程對比Fig.13 Comparison of landslide processes

4 結論

1)采用Geostudio與PFC方法耦合計算降雨工況下邊坡變形穩(wěn)定問題，可同時克服有限元方法難以計算邊坡大變形及離散元方法滲流計算效率較低的缺陷，通過該耦合方法可分析降雨工況下邊坡表面沉降、坡面沖刷、土體崩塌等現(xiàn)象。

2)對于單一介質均質邊坡，降雨入滲規(guī)律與土體的所在位置無關，僅與其所在深度相關，沿降雨入滲方向，會形成一定厚度的飽和帶和非飽和帶。

3)大雨工況下，邊坡變形主要為土體沉降，變形最大區(qū)域位于邊坡頂部；暴雨工況下，由于坡面徑流的沖刷作用，坡腳處首先發(fā)生破壞，隨后誘發(fā)了后部土體大范圍滑動，發(fā)生漸進式牽引式滑坡。