基于ABAQUS強度折減法分析庫水位下降對邊坡穩定性影響

郭 飛 付調金 阮榮乾 管宏飛

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

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

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

1 樹坪滑坡概況

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

圖1 樹坪滑坡全貌

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

圖2 樹坪滑坡工程地質剖面圖

滑體2為滑坡堆積物（Qdel），主要為碎塊石粘土層，碎塊石成分主要為泥質粉砂巖、泥灰巖和灰巖.土為壤土、粉質粘土，填充于碎塊石中，結構不均勻，稍密，滑體堆積物結構松散，透水性較好.

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

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

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

2.1 滑坡飽和與非飽和滲流場與應力場耦合模型［4]

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

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

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

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

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

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

2.2 計算模型確定

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

圖3 樹坪滑坡計算模型

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

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

2.3 計算工況

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

圖4 三峽水庫運行水位調度圖

2.4 計算參數選取

2.4.1 非飽和土滲透性參數

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

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

表1 樹坪滑坡1號滑體Sr－Uw對應關系 （孔壓單位：kPa）

表2 樹坪滑坡1號滑體Sr－Ks對應關系

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

表3 樹坪滑坡2號滑體Sr－Uw對應關系 （孔壓單位：kPa）

表4 樹坪滑坡2號滑體Sr－Ks對應關系

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

表5 樹坪滑坡滑帶Sr－Uw對應關系 （孔壓單位：kPa）

表6 樹坪滑坡滑帶Sr－Ks對應關系

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

2.4.2 材料力學參數

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

表7 樹坪滑坡巖土體物理學力學參數

3 ABAQUS強度折減法原理

3.1 強度折減法原理［5]

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

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

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

3.2 失穩判據

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

2）以特征部位的位移拐點作為評價標準.根據計算域內某一部位處結點的位移與折減系數之間關系密切曲線的變化特征確定失穩狀態，當某處一結點位移突然迅速增大，則認為邊坡發生失穩［8].

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

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

4 結果分析

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

圖5 樹坪滑坡不同時刻滑體內塑性區分布圖

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

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

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

5 結 論

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

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

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

［2]黃潤秋，戚國慶.非飽和滲流基質吸力對邊坡穩定性的影響［J].工程地質學報，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流固耦合理論的庫岸滑坡穩定性分析［D].濟南：山東大學，2005.

［5]費 康，張建偉.ABAQUS在巖土工程中的應用［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]呂擎峰.土坡穩定分析方法研究［D].南京：河海大學，2005.

［8]宋二祥.土工結構安全系數的有限元計算［J].巖土工程學，1997，19（2）：1－7.

［9]欒茂田，武亞軍，年廷凱.強度折減有限元方法中邊坡失穩的區判據及其應用［J].防災減災工程學報，2003，23（3）：1－8.