李子隆，許 模，趙 勇，王子忠，郭 健，李明月

(1.成都理工大學(xué)地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點實驗室，四川成都610059；2.四川省水利水電勘測設(shè)計研究院，四川成都610072)

0 引 言

隨著人類對水資源需求的不斷增加，水利水電工程的數(shù)量也逐漸增長，庫岸邊坡的安全問題更不容忽視。庫岸邊坡的失穩(wěn)會帶來嚴(yán)重的災(zāi)難，大量的巖土體進(jìn)入水庫產(chǎn)生涌浪(如Vaiont水庫)，巖土體堆積造成堵江并形成堰塞湖，對人類的生存環(huán)境造成嚴(yán)重威脅。目前，對邊坡穩(wěn)定性的研究方法層出不窮，數(shù)值模擬軟件各式各樣，但對滑坡失穩(wěn)后造成堵江的研究通常都是采用經(jīng)驗法[1-2]，相關(guān)研究軟件也相對較少，最常用的方法是用FLAC3D分析邊坡的穩(wěn)定性及變形特征，再用UDEC模擬邊坡的破壞特征并預(yù)測堰塞壩高度[3- 4]。徐文杰對肖家河滑坡穩(wěn)定性做了分析，并用ABAQUS/Explicit動力有限元分析滑坡在地震工況下的堵江過程[5]。有限元模擬的塊體不能破裂，模擬效果并不明顯。學(xué)者們對于土質(zhì)邊坡的模擬大多采用離散元顆粒流程序PFC。宋朋燃考慮了降雨的因素，用PFC3D模擬了黃土滑坡坡面沖刷破壞過程[6]；馬建全用PFC2D模擬了黑方臺滑坡各破壞階段的變形特征；周健建立了粘性土坡和砂性土坡的理想模型，分析了微觀參數(shù)對斜坡變形的影響，隨著參數(shù)的增加，變形由塑性轉(zhuǎn)變?yōu)榇嘈訹7]；梁承洋用PFC2D模擬了川藏交通廊道冰磧物滑坡堵江堆積物的形態(tài)，同時研究了單因素變量對堵江的影響[8]。

地下水是庫岸邊坡失穩(wěn)的最重要因素，很多欠穩(wěn)定的斜坡在地下水位升高后覆蓋層失穩(wěn)，并造成堵江，但對于考慮地下水位的邊坡堵江研究并不多。為此，本文采用PFC2D顆粒流程序，建立有覆蓋層的庫岸邊坡模型，采用降低微觀參數(shù)的方法模擬不同的地下水位，并研究其破壞過程和堵江壩高預(yù)測，分析庫岸邊坡覆蓋層飽和土體比例與堵江形成堰塞壩高度之間的關(guān)系。

1 顆粒流PFC理論

1.1 基本假設(shè)

PFC2D是基于離散單元法的顆粒流程序軟件，通過顆粒的位移和接觸構(gòu)成整個系統(tǒng)的計算單元。模型中主要包括球單元(ball)和墻單元(wall)，以及顆粒與顆粒接觸和顆粒與墻接觸。在PFC2D中，做了如下假設(shè)[9-12]：

(1)顆粒被視為剛體，不可變形。

(2)顆粒之間的接觸是點接觸，范圍小。

(3)顆粒接觸為柔性接觸，允許一定的重疊，但重疊量相對于單元尺寸足夠小。

(4)顆粒間可以存在粘結(jié)，且遵循力-位移定律。

軟件采用顯式中心差分法迭代計算顆粒的速度和加速度。每個計算步長內(nèi)顆粒只與接觸的顆粒相互作用，所有的顆粒遵循牛頓第二定律，墻不遵循牛頓第二定律。每次計算后更新球的位置，通過運動定律和力-位移定律反復(fù)迭代計算，且2種定律的計算同時進(jìn)行。計算過程見圖1。

圖1 PFC2D循環(huán)計算

1.2 運動方程

力-位移方程描述了顆粒間的力和位移關(guān)系，顆粒間法向力-位移方程為

ΔFn=knun

式中，F(xiàn)n為法向接觸力；kn為法向接觸剛度；un為相對位移。

顆粒間切向力-位移方程為

ΔFs=ksus

式中，F(xiàn)s為切向接觸力；ks為切向接觸剛度；us為相對位移。

顆粒的平動和轉(zhuǎn)動是通過力和力矩計算平動加速度和轉(zhuǎn)動加速度，再根據(jù)加速度計算顆粒在時間Δt內(nèi)的平動位移量和轉(zhuǎn)動位移量，以x方向為例，時間t0顆粒的平動加速度和轉(zhuǎn)動加速度分布為

在時間t2=t0+Δt時，顆粒的位移為

式中，ux為x方向的位移。

1.3 連接模型

PFC2D顆粒連接模型包括接觸點連接模型和接觸平行連接模型。點連接模型只能傳遞力，用法向/切向連接強度定義，但不能傳遞力矩。當(dāng)法/切向力大于法/切向連接強度時，連接破壞。而平行連接模型可以視為以接觸點為中心的一系列均勻分布且相互平行的彈簧，可以傳遞力矩，用平行連接法向/切向強度定義。當(dāng)法/切向應(yīng)力大于法/切向連接強度時，連接破壞。對于庫岸邊坡覆蓋層堆積物，選擇平行連接模型更加吻合。

1.4 參數(shù)標(biāo)定

圖2 宏觀參數(shù)與微觀參數(shù)關(guān)系

法向接觸剛度kn/N·m-1切向接觸剛度ks/N·m-1摩擦系數(shù)f非飽和狀態(tài)平行法向連接強度/Pa·m-1平行切向連接強度/Pa·m-1密度/kg·m-3飽和狀態(tài)平行法向連接強度/Pa·m-1平行切向連接強度/Pa·m-1密度/kg·m-36×1086×1080.69×1059×10519606.5×1056.5×1052100

圖5 各水位模擬結(jié)果

2 模型建立及結(jié)果分析

2.1 模型建立

由于水庫的運行，在庫岸邊坡內(nèi)會產(chǎn)生不同高度的地下水位，對覆蓋層穩(wěn)定性有較大的影響。建立的較理想的邊坡模型見圖3。覆蓋層水平長200 m，縱向高差50 m，基覆面傾角約30°，河谷寬40 m，深9 m。假設(shè)在天然水位下覆蓋層處于基本穩(wěn)定狀態(tài)，水庫在運行工況下分別產(chǎn)生了高于河谷5、10、15、20、25 m和30 m的不同高度的地下水位。覆蓋層的體積約為2 000 m3，天然水位下的覆蓋層PFC模型見圖4。本次計算模型共2 196個顆粒，粒徑0.25～0.75 m。模型選用平行連接結(jié)構(gòu)，在對微觀參數(shù)進(jìn)行調(diào)整之后，使模型處于欠穩(wěn)定的臨界狀態(tài)。微觀參數(shù)取值見表1。在不同水位高度下，飽和土體的比例見表2。各水位模擬結(jié)果見圖5。

圖3 理想的庫岸邊坡覆蓋層(單位：m)

圖4 天然水位下的覆蓋層PFC模型

水位高度/m飽和土比例51∶20101∶5151∶3201∶2253∶5304∶5

2.2 結(jié)果分析

從圖5可知，總體上，覆蓋層破壞主要是由于底部土體在地下水的作用下力學(xué)參數(shù)降低，在河流側(cè)蝕作用下發(fā)生局部破壞，進(jìn)一步導(dǎo)致覆蓋層臨空狀態(tài)的改變，再接著引起新的臨空面產(chǎn)生，后續(xù)的破壞將繼續(xù)進(jìn)行。在水位為5、10 m和20 m時，覆蓋層的破壞模式是塌岸，造成臨空狀態(tài)改變，但中～后段的土體并沒有因此而產(chǎn)生較大的破壞。從趨勢上看，水位越高，底部土體塌岸的規(guī)模越大，中～后段的變形量越大。在水位為25 m時，前～中段的土體發(fā)生了失穩(wěn)，雖然后段的土體沒有大量的位移，但也產(chǎn)生了新的較大的拉裂。在水位30 m時，覆蓋層絕大部分被地下水浸沒，整個覆蓋層內(nèi)部的力學(xué)強度降低，造成整體破壞。

斜坡失穩(wěn)后造成的堵江高度見圖6。從圖6可知，在水位高度為15 m時，模擬的堵江高度是本次研究當(dāng)中最大的，其破壞過程見圖7。圖中，T為運算步數(shù)。從圖7可知，在初始時期，覆蓋層的破壞主要是地下水位以下部分；中期階段，飽和土體的下滑對上部非飽和土體產(chǎn)生一個拖拽力，主要發(fā)生在地下水位線附近，帶動非飽和土體一起下滑，由于不斷產(chǎn)生新的臨空面，后面的土體不能達(dá)到平衡狀態(tài)，因此也跟著下滑，最終導(dǎo)致斜坡的整體失穩(wěn)。

圖6 地下水位高度與堵江壩高關(guān)系

圖7 水位為15 m時斜坡破壞過程

用同樣的方法分析其他地下水位高度時的破壞過程。在5 m和10 m時，由于坡腳的飽和土體下滑，產(chǎn)生新的臨空面，但中～后段大部分土體未受到影響，其基本結(jié)構(gòu)并未遭到改變，新形態(tài)的斜坡尚能達(dá)到平衡狀態(tài)。在20 m和25 m時，斜坡失穩(wěn)的主要部位是飽和土體，在飽和土體下滑后，非飽和土體位移并不明顯，且在飽和土體下滑部位延伸至地下水位線時，并未造成飽和土體和非飽和土體一起大量下滑(同15 m時的情況)。其原因為上部非飽和土體并沒有15 m時多，飽和土體能夠支撐住非飽和土體，不使其進(jìn)一步下滑。而30 m時，由于飽和土體下滑較多，造成了后段斜坡的不穩(wěn)定。通過對比各地下水位高度時，斜坡覆蓋層破壞后的堵江高度可知，整體上，地下水位越高，堵江堰塞壩高度越高，但在飽和土體的比例占1/3時，壩高比其相鄰水位高出很多，即存在一個最危險水位。

3 結(jié) 語

本文基于PFC2D顆粒流程序，建立有覆蓋層的庫岸邊坡模型，得出以下結(jié)論：

(1)整體上，地下水位越高，斜坡覆蓋層失穩(wěn)后堆積形成的堰塞壩高度越高。但在飽和土體的比例占1/3時，壩高比其相鄰水位高出很多，即存在一個最危險水位。

(2)水位低于危險水位時，斜坡覆蓋層主要破壞模式是坡腳塌岸；水位高于危險水位區(qū)域時，主要是飽和土體發(fā)生失穩(wěn)，非飽和土體部分失穩(wěn)并堆積在水位線附近。

雖然PFC2D能夠較好地模擬土質(zhì)斜坡堵江，但土體的微觀參數(shù)取值較復(fù)雜，在有試驗數(shù)據(jù)的條件下，可用雙軸試驗來標(biāo)定參數(shù)；在沒有試驗數(shù)據(jù)時，則需類比和多次試模擬以確定微觀參數(shù)。此外，靜態(tài)水位可用PFC2D模擬，只能用于水庫調(diào)動之后的水位穩(wěn)定的情況，而對庫水位升降造成的地下水滲流模擬難度較大。

參考文獻(xiàn)：

[1] 孫蔚, 董云. 堵江滑坡壩的潰壩方式模擬與環(huán)境效應(yīng)分析[J].自然災(zāi)害學(xué)報, 2008, 17(4)： 72- 77．

[2] 馬建全. 黑方臺灌區(qū)臺緣黃土滑坡穩(wěn)定性研究[D]. 長春： 吉林大學(xué), 2012．

[3] 原俊紅. 白龍江中游滑坡堵江問題研究[D]. 蘭州： 蘭州大學(xué), 2007．

[4] 羅剛. 唐家山高速短程滑坡堵江及潰壩機制研究[D]. 成都： 西南交通大學(xué), 2012．

[5] 徐文杰, 陳祖煜, 何秉順, 等. 肖家橋滑坡堵江機制及災(zāi)害鏈效應(yīng)研究[J]. 巖石力學(xué)與工程學(xué)報, 2010, 29(5)： 933- 942．

[6] 宋朋燃. 黃土邊坡沖刷破壞特征及數(shù)值模擬[D]. 長春： 吉林大學(xué), 2013．

[7] 周健, 王家全, 曾遠(yuǎn), 等. 土坡穩(wěn)定分析的顆粒流模擬[J]. 巖土力學(xué), 2009, 30(1)： 86- 90．

[8] 梁承洋. 川藏交通廊道冰磧物滑坡堵江風(fēng)險及線路對策研究[D]. 成都： 西南交通大學(xué), 2015．

[9] 周健, 池永, 池毓蔚, 等. 顆粒流方法及PFC2D程序[J]. 巖土力學(xué), 2000, 21(3)： 1- 4．

[10] CUNDALL P A, STRACK O D L. A discrete numerical mode for granular assemblies[J]. Géotechnique, 1979, 29(1): 47- 65．

[11] 彭博, 陳玉明, 袁利偉. 基于 PFC2 D馬鞍山鐵礦2號露天采場邊坡穩(wěn)定性分析[J]. 黃金, 2014(12)： 37- 40．

[12] 陳俊, 張東, 黃曉明. 離散元顆粒流軟件(PFC)在道路工程中的應(yīng)用[M]. 北京： 人民交通出版社股份有限公司, 2015．