郭 飛 付調(diào)金 阮榮乾 管宏飛

（1.三峽大學 三峽庫區(qū)地質(zhì)災(zāi)害教育部重點實驗室，湖北 宜昌 443002；2.長江三峽勘測研究院有限公司，武漢 430074；3.水利部海委漳衛(wèi)南運河樂陵河務(wù)局，山東 德州 253600）

水是誘發(fā)滑坡的主要因素.據(jù)不完全統(tǒng)計，三峽庫區(qū)在175m水位范圍內(nèi)共有大小滑坡2 000余個.在日本，大約60%的水庫滑坡發(fā)生在庫水位驟降時期，40%發(fā)生在水位上升時期，包括蓄水初期［1].庫水位驟然下降時，由于坡體內(nèi)地下水位下降相對滯后，導致坡體內(nèi)產(chǎn)生超孔隙水壓力，對滑坡的穩(wěn)定性會產(chǎn)生不利影響，因此，三峽庫區(qū)水位下降對滑坡穩(wěn)定性的影響是當前滑坡研究領(lǐng)域十分重要的研究課題.

目前對三峽庫區(qū)滑坡穩(wěn)定性評價和治理都是建立在飽和土假設(shè)基礎(chǔ)上，對暴雨、庫水位下降產(chǎn)生滑坡的機理研究較少，已有的工程設(shè)計均采用對暫態(tài)飽和區(qū)及暫態(tài)水壓力進行假定的方法［2].同時，傳統(tǒng)的極限平衡法難以有效地考慮水與巖土體的相互作用機理以及滲流場和應(yīng)力場的耦合效應(yīng).基于飽和與非飽和的流固耦合理論的數(shù)值方法雖然能解決這些問題，但它是根據(jù)邊坡的位移場、應(yīng)力場、塑性區(qū)來間接地評價，不能直接給出安全系數(shù)以直觀評價邊坡穩(wěn)定性.本文通過ABAQUS軟件，考慮飽和與非飽和流固耦合效應(yīng)的強度折減分析，可以直接獲得一個安全系數(shù)，此方法不僅保持了有限元在模擬復(fù)雜問題上的優(yōu)點，而且概念明確，結(jié)果直觀，對實際邊坡治理工程具有一定的參考意義.

1 樹坪滑坡概況

樹坪滑坡［3]位于湖北省秭歸縣沙鎮(zhèn)溪鎮(zhèn)長江干流南岸斜坡地段的古滑坡.上距沙鎮(zhèn)溪鎮(zhèn)千將坪滑坡約3km處，下距三峽工程三斗坪壩址47km.樹坪滑坡屬多期性巨型滑坡，它是由兩個區(qū)塊（見圖1中區(qū)塊1和區(qū)塊2）組成，物質(zhì)組成較復(fù)雜（見圖2），樹坪滑坡區(qū)塊1大致可以分為滑體、滑帶、滑床3部分，其中滑體有兩種組成物質(zhì).

圖1 樹坪滑坡全貌

滑體1為坡積物（Qdl），主要為粉質(zhì)粘土層夾碎塊石，碎塊石多呈次棱角狀，以粉砂巖、泥質(zhì)粉砂巖為主，結(jié)構(gòu)松散～稍密，透水性較差.

圖2 樹坪滑坡工程地質(zhì)剖面圖

滑體2為滑坡堆積物（Qdel），主要為碎塊石粘土層，碎塊石成分主要為泥質(zhì)粉砂巖、泥灰?guī)r和灰?guī)r.土為壤土、粉質(zhì)粘土，填充于碎塊石中，結(jié)構(gòu)不均勻，稍密，滑體堆積物結(jié)構(gòu)松散，透水性較好.

滑帶為紫紅色角礫石土層，較濕，結(jié)構(gòu)緊密，土可塑.碎石呈次棱角狀～次園狀，碎石上可見擦痕，土層中可見明顯揉皺、光滑鏡面.滑帶埋深較大，層厚一般在10～20cm.

自2003年6月三峽水庫第一次蓄水至高程135 m后，庫區(qū)發(fā)生大量的滑坡復(fù)活及失穩(wěn)現(xiàn)象.2004年1月，樹坪滑坡出現(xiàn)多處變形，宏觀變形以地面裂縫和房屋開裂為主.根據(jù)地質(zhì)勘察結(jié)果，坡體側(cè)緣剪切帶內(nèi)出現(xiàn)新鮮錯位，江面以下剪出口處，特別是在降雨時節(jié)，江面出現(xiàn)大面積渾水現(xiàn)象.但根據(jù)自動位移計的監(jiān)測資料表明，目前樹坪滑坡處于較穩(wěn)定狀態(tài).

2 庫水下降飽和與非飽和流固耦合數(shù)值模擬

2.1 滑坡飽和與非飽和滲流場與應(yīng)力場耦合模型［4]

由應(yīng)力場平衡方程（1）和非飽和滲流場連續(xù)方程（2）共同組成了流固耦合問題的控制方程組：

其中，方程（1）中IN是內(nèi)力矩陣，IP是外力矩陣，此方程由虛功原理推導簡化而來.

流固耦合數(shù)學模型的求解，還需要有相應(yīng)的定解條件，即確定模型的邊界條件和初始條件才能求解.本文研究庫水下降對滑坡的影響作用，其邊界條件需滿足方程（3）和方程（4）：

其中，方程（3）為水頭邊界條件，方程（4）為初始孔隙水壓力的分布函數(shù).

該方程組還需同時滿足以下初始條件：位移邊界條件，即U＝u（x，y，z）；初始時刻（t＝0），位移或質(zhì)點速度的初始值.

將流固耦合控制方程組（1）、（2）與在庫水下降時相應(yīng)的滲流場與應(yīng)力場邊界條件和初始條件相結(jié)合，就可得到解決非飽和流固耦合問題的控制方程，然后便可求解方程組.

2.2 計算模型確定

根據(jù)樹坪滑坡的野外地質(zhì)勘察和鉆孔資料，選擇計算剖面如圖3所示，該剖面為樹坪滑坡1號塊體的主滑方向，沿該方向上布置有日本坂田電機株式會社生產(chǎn)的自動位移計，計算模型尺寸為：X方向上最大長度為946m（0＜X＜946m），模型后緣Y方向上最大高程為468m（0＜Y＜468m），前緣Y方向上45m.采用ABAQUS數(shù)值計算軟件進行工況模擬，計算域包含滑體、滑帶和基巖，整個計算域剖分網(wǎng)絡(luò)單元12 390個，結(jié)點共25 330個.

圖3 樹坪滑坡計算模型

約束邊界條件：左右邊界約束水平位移，下邊界約束水平和豎直位移，上邊界自由.

滲流邊界條件：右側(cè)邊界庫水位以下為水頭邊界，庫水位以上為零流量邊界；左側(cè)邊界為水頭邊界.初始地下水位線根據(jù)地質(zhì)鉆孔資料得到的實際地下水位線，并通過兩側(cè)水頭計算得出.

2.3 計算工況

根據(jù)三峽水庫運行曲線圖（見圖4），每年5月至9月為汛期，庫水位保持145m低水位運行；9月底至11月初，汛期過后為滿足發(fā)電需要，水位從145m大幅上升至175m；11月至次年1月為枯水期，水位維持在175m高水位運行；1月至5月為騰出防洪庫容，水位緩慢降落至145m.目前，較多學者研究庫水位上升及驟降對庫區(qū)邊坡的影響，未充分考慮到從1～5月，水位緩慢降落情況，從水位調(diào)度圖中可以看出三峽庫區(qū)在此時間段內(nèi)庫水位下降平均速度約為0.21m／d，而大多數(shù)學者考慮水位驟降的影響，對于揭示古滑坡在庫水位變動作用下的復(fù)活機理有一定意義，但不能反映水庫正常運行下對古滑坡的影響.本文考慮水位正常調(diào)度情況下，庫水位下降（下降速度約為0.21m／d）對樹坪滑坡的影響.

圖4 三峽水庫運行水位調(diào)度圖

2.4 計算參數(shù)選取

2.4.1 非飽和土滲透性參數(shù)

非飽和滲流分析涉及的參數(shù)除了飽和滲透系數(shù)外，還需要材料的初始含水率、土水特征曲線及非飽和滲透性函數(shù).由于樹坪滑坡非飽和土性狀還未進行實驗研究，則其滑體和滑帶的非飽和參數(shù)采用工程類比法，同時根據(jù)地質(zhì)勘察資料確定樹坪滑坡3種材料的滲透性參數(shù)分別如下：

1）1號滑體：滲透系數(shù)取為2.8×10－5m／s，即2.419m／d.飽和度（Sr）與孔隙水壓力（Uw）及滲透系數(shù)折減系數(shù)（Ks）關(guān)系見表1～2.

表1 樹坪滑坡1號滑體Sr－Uw對應(yīng)關(guān)系 （孔壓單位：kPa）

表2 樹坪滑坡1號滑體Sr－Ks對應(yīng)關(guān)系

2）2號滑體：滲透系數(shù)取為6.0×10－5m／s，即5.184m／d.飽和度（Sr）與孔隙水壓力（Uw）及滲透系數(shù)折減系數(shù)（Ks）關(guān)系見表3～4.

表3 樹坪滑坡2號滑體Sr－Uw對應(yīng)關(guān)系 （孔壓單位：kPa）

表4 樹坪滑坡2號滑體Sr－Ks對應(yīng)關(guān)系

3）滑帶：滲透系數(shù)取為5.0×10－7m／s，即0.043 m／d.飽和度（Sr）與孔隙水壓力（Uw）及滲透系數(shù)折減系數(shù)（Ks）關(guān)系見表5～6.

表5 樹坪滑坡滑帶Sr－Uw對應(yīng)關(guān)系 （孔壓單位：kPa）

表6 樹坪滑坡滑帶Sr－Ks對應(yīng)關(guān)系

4）滑床：滲透系數(shù)取為2.0×10－10m／s，即1.728×10－5m／d.因不考慮基巖的非飽和特性，所以，計算參數(shù)不含飽和度（Sr）與孔隙水壓力（Uw）及滲透系數(shù)折減系數(shù)（Ks）關(guān)系.

2.4.2 材料力學參數(shù)

應(yīng)力應(yīng)變分析采用 Mohr－Coulomb非飽和土本構(gòu)模型，計算參數(shù)主要有巖土體變形模量（E）、泊松比（μ）、容重（γ）、凝聚力（c）、內(nèi)摩擦角（φ），主要由地質(zhì)勘察資料所建議.滑坡計算參數(shù)見表7.

表7 樹坪滑坡巖土體物理學力學參數(shù)

3 ABAQUS強度折減法原理

3.1 強度折減法原理［5]

強度折減法最早由Zienkiewicz等提出，后被許多學者廣泛采用.他們提出了一個抗剪強度折減系數(shù)（SSRF：Shear Strength Reduction Factor）的概念，其定義為：在外載荷作用保持不變的情況下，外載荷所產(chǎn)生的實際剪應(yīng)力與抵御外載荷所發(fā)揮的最低抗剪強度是按照實際強度指標折減后所確定的、實際中得以發(fā)揮的抗剪強度相等.當假定邊坡內(nèi)所有土體抗剪強度的發(fā)揮程度相同時，這種抗剪強度折減系數(shù)相當于傳統(tǒng)意義上的邊坡整體穩(wěn)定安全系數(shù)Fs，又稱為強度儲備安全系數(shù)，與極限平衡法中所給出的穩(wěn)定安全系數(shù)在概念上是一致的.

折減后的抗剪強度參數(shù)可分別表示為

式中，c和φ是土體所能夠提供的抗剪強度參數(shù)；cm和φm是維持平衡所需要的或土體實際發(fā)揮的抗剪強度；Fr是強度折減系數(shù).

3.2 失穩(wěn)判據(jù)

1）以數(shù)值計算收斂與否作為評價標準.即根據(jù)有限元解的收斂性確定失穩(wěn)狀態(tài).該失穩(wěn)判據(jù)與一定的計算方法相關(guān).因此，以數(shù)值計算收斂與否為失穩(wěn)判據(jù)是不合理的，適用性差［6－7].

2）以特征部位的位移拐點作為評價標準.根據(jù)計算域內(nèi)某一部位處結(jié)點的位移與折減系數(shù)之間關(guān)系密切曲線的變化特征確定失穩(wěn)狀態(tài)，當某處一結(jié)點位移突然迅速增大，則認為邊坡發(fā)生失穩(wěn)［8].

3）以是否形成連續(xù)的貫通區(qū)作為評價標準.理論上，邊坡的變形過程總是伴隨著一些物理量的出現(xiàn)和發(fā)展，如塑性區(qū)、塑性應(yīng)變、廣義剪應(yīng)變和應(yīng)力水平等，當這些物理量達到一定的值時，邊坡失穩(wěn)［9].

本文采用強度折減法分析庫水下降對樹坪滑坡影響時，主要采用第2）和3）兩種判據(jù)，分析在庫水下降過程中，滑坡的位移場、塑性區(qū)及穩(wěn)定系數(shù)，具有較為明確的物理意義.

4 結(jié)果分析

通過ABAQUS軟件，采用彈塑性有限元強度折減法對樹坪滑坡進行飽和－非飽和流固耦合計算，其計算結(jié)果如圖5～6所示.

圖5 樹坪滑坡不同時刻滑體內(nèi)塑性區(qū)分布圖

圖6 樹坪滑坡位移增量云圖（單位：m）

由圖5可以看出，庫水由175m以0.21m／d的速度下降至145m過程中，塑性應(yīng)變主要發(fā)生在滑帶的前緣，滑帶的后緣和中部基本上沒有產(chǎn)生塑性應(yīng)變.最大的塑性應(yīng)變發(fā)生在滑帶的前緣，對應(yīng)著滑坡的剪出口.而且，隨著庫水不斷下降，塑性應(yīng)變有擴大的趨勢.根據(jù)地質(zhì)資料分析，在庫水作用下，樹坪滑坡變形破壞方式主要是牽引式解體，慢速滑移.從塑性區(qū)的分布來看，計算結(jié)果基本上和地質(zhì)上的認識和判斷是較符合的.樹坪滑坡之所以未產(chǎn)生整體性滑移，主要可能由于滑帶中部提供強有力的抗滑力，如果滑帶中部產(chǎn)生較大的塑性變形，則會導致樹坪滑坡滑帶塑性區(qū)全部貫通，進而導致其整體性失穩(wěn).

從位移云圖分析得出，樹坪滑坡的最大位移發(fā)生在滑體1和滑體2前緣位置.從滑體前緣至后緣，位移逐漸減小，后緣和中部基本無位移，與塑性區(qū)的分布吻合，此時滑坡的位移主要由前緣帶動.當采用特征部位的位移拐點判據(jù)得到的安全系數(shù)為1.112，此時滑坡處于較穩(wěn)定狀態(tài).

5 結(jié) 論

1）通過數(shù)值計算，在庫水以0.21m／d的速度由175m降至145m過程中，樹坪滑體1和滑體2前緣出現(xiàn)較大位移，滑體后緣和中部基本上未產(chǎn)生位移，計算結(jié)果表明，樹坪滑坡由于前緣的牽引，慢速滑移.這與地質(zhì)資料較符合.

2）通過數(shù)值模擬，三峽水庫在正常泄水的過程中，樹坪滑坡塑性區(qū)和位移出現(xiàn)的位置，表明在該工況下，樹坪滑坡的穩(wěn)定性受到了一定的影響，引起滑體前緣變形，但計算穩(wěn)定性系數(shù)為1.112，與監(jiān)測資料較吻合.該計算結(jié)果為三峽水庫的調(diào)度及災(zāi)害防治提供一定的參考.

［1]朱冬林，任光明，聶德新，等.庫水位變化下對水庫滑坡穩(wěn)定性影響的預(yù)測［J].水文地質(zhì)工程地質(zhì)，2002，3（2）：6－9.

［2]黃潤秋，戚國慶.非飽和滲流基質(zhì)吸力對邊坡穩(wěn)定性的影響［J].工程地質(zhì)學報，2002，10（4）：343－348.

［3]Wang Fawu，Zhang Yemin，et al.Deformation Features of Shuping Landslide Caused by Water Level Changes［J].Chinese Journal of Rock Mechanics and Engineering，2007，26（3）：509－517.

［4]張 欣.基于ABAQUS流固耦合理論的庫岸滑坡穩(wěn)定性分析［D].濟南：山東大學，2005.

［5]費 康，張建偉.ABAQUS在巖土工程中的應(yīng)用［M].北京：中國水利水電出版社，2010.

［6]Dawson E M，Roth W H，Drescher A.Slope Stability Analysis by Strength Reduction［J].Geotechnique，1999，49（3）：835－840.

［7]呂擎峰.土坡穩(wěn)定分析方法研究［D].南京：河海大學，2005.

［8]宋二祥.土工結(jié)構(gòu)安全系數(shù)的有限元計算［J].巖土工程學，1997，19（2）：1－7.

［9]欒茂田，武亞軍，年廷凱.強度折減有限元方法中邊坡失穩(wěn)的區(qū)判據(jù)及其應(yīng)用［J].防災(zāi)減災(zāi)工程學報，2003，23（3）：1－8.