基于孔壓-應力耦合下考慮滲透系數(shù)各向異性對庫水位驟降邊坡的影響

戚海棠，郁舒陽，任旭華，張繼勛

(河海大學水利水電學院，江蘇 南京 210098)

庫水位的上升導致滑坡土體內部孔隙水壓增大，土體的強度參數(shù)及有效應力降低[1-3]；此外，庫水位的下降導致滑坡內部孔壓變化滯后庫水位，對周圍土體骨架產生向臨空面及向下的滲流動水壓力(動態(tài)擴張力)，牽引帶動坡面土體[4-6]，因此，對庫水位變動下的滑坡失穩(wěn)規(guī)律及影響因素的把握是正確認識庫水位滲流作用機理及防治滑坡災害的關鍵所在。

國內外學者針對庫水位波動下的滑坡滲流力學響應展開了大量研究，其中，IQBAL J等[7]分析研究滑坡與周圍環(huán)境因素的關系，LUO F Y等[8]研究發(fā)現(xiàn)隨著坡度的增加、庫水位下降下的邊坡滑移面逐漸變深，HU X L等[9]利用非飽和滲流理論及有限元原理對黃土坡濱江滑坡穩(wěn)定性進行了數(shù)值模擬分析，JIAN W X等[10]對不同降雨強度及庫水位驟降速率下的千江坪滑坡的滲流穩(wěn)定性進行了數(shù)值模擬，但是這些研究未考慮邊坡土體的滲流各向異性。宋云奇等[11]研究表明黏土等土質因其絮狀微觀結構導致滲透系數(shù)存在各向異性，同時土質的滲透系數(shù)各向異性受其干密度及凍融循環(huán)等因素的影響巨大[12]；YEH H F[13]考慮土質邊坡的滲透系數(shù)各向異性程度，對其滲流特性及局部安全系數(shù)進行數(shù)值模擬，但忽略了各向異性方向的影響。

由2013年6月至2015年12月三峽庫區(qū)的水位波動及滑坡數(shù)量統(tǒng)計示意圖[14-15]可知，三峽庫區(qū)庫水位變動是導致三峽庫區(qū)滑坡的重要誘因。針對以往研究的不足，本文以三峽庫區(qū)蔡坡堆積體為例[16]，對砂土與黏土2種滑坡土質，考慮到庫水位下降下不同滲透系數(shù)各向異性程度kr及方向α情況的滲流、變形及穩(wěn)定性進行有限元分析，從而為全面認識邊坡滲流各向異性規(guī)律及對滑坡災害的防治提供了一定的參考。

1 非飽和孔壓-應力耦合計算模型

非飽和滲流各向異性孔壓-應力耦合及穩(wěn)定性計算利用Geo-studio中的SIGMA/W與SLOPE/W模塊進行耦合分析，SIGMA/W對有限元網(wǎng)格的每一個節(jié)點創(chuàng)建3個方程，其中二個是平衡方程(位移)，另一個是滲流連續(xù)方程(孔隙水壓力)，同時，求解3個方程能夠認識位移跟孔隙水壓力的變化。SLOPE/W可以根據(jù)SIGMA/W計算得出的孔壓進行非飽和安全系數(shù)的計算。

1.1 非飽和土滲流有限元計算理論

根據(jù)飽和及非飽和達西定律[17]推導出滲流計算有限元公式中的滲流控制方程，其二維表達式為

(1)

式(1)中x、y分別為x、y方向的位置坐標，kx為x方向的滲透系數(shù)，ky為y方向的滲透系數(shù)，H為總水頭，Q為施加的邊界流量，t為時間，mw為儲水曲線的斜率，γw為水的容重。

1.2 滲流-應力耦合計算理論

對于二維問題，非飽和土介質的增量應變-應力關系可以表達成為：

{Δσ}=[D]{Δε}-[D]{mH}(ua-uw)+{Δua},

(2)

式(2)中ε為法向應變，γ為工程剪應變，σ為法向應力，τ為剪切應力，ua為孔隙氣壓力，uw為孔隙水壓力，[D]為排水本構矩陣，{mH}可以表達為：

(3)

式(3)中H′為與基質吸力(ua-uw)有關的函數(shù)。

1.3 非飽和安全系數(shù)

邊坡安全計算系數(shù)采用基于極限平衡理論的Morgenstern-Price法，此法嚴格滿足力平衡及力矩平衡，計算精度較高。

2 計算模型及邊界條件

計算模型選擇三峽庫區(qū)蔡坡堆積體[16]為例，蔡坡堆積體位于鶴峰縣燕子鄉(xiāng)境內，距壩址上游約為1.32 km，其平面形態(tài)近似矩形，順岸坡展布，地形坡角15°～40°，前緣高程約160 m，后緣高程220 m，總體東高西低。滑坡厚度為2～4 m，總面積約1.2萬m2，體積約3.6萬m3。為了減小邊界條件對計算結果的影響，將坡腳和坡頂?shù)姆秶M行延伸，整個模型共劃分為1 084個節(jié)點、1 001個單元,計算模型如圖1所示。

圖1 計算模型及邊界示意圖

為了綜合研究邊坡土體各向異性對邊坡不同位置的滲流、變形特性的影響，設置上部(x=100 m)、中部(x=200 m)、下部(x=300 m)3個監(jiān)測面反映邊坡不同位置的滲流、變形特性，監(jiān)測面貫穿滑坡體。初始條件為AB 190 m、DF 175 m水頭計算所得的穩(wěn)定滲流場；邊界條件設置如下：AB為190 m水頭邊界，DEF為175～145 m庫水位變動邊界，AF、CD為不透水邊界。

3 材料參數(shù)及計算工況

3.1 非飽和計算參數(shù)

土-水特征曲線(SWCC)采用Fredlund & Xing模型[18]進行估算，該方法的控制方程見文獻[18]。

邊坡材料選取典型的砂土(滲透系數(shù)較大)及黏土(滲透系數(shù)較小)[19]，進行邊坡非飽和滲流分析時的計算參數(shù)設置[20]如表1所示，同時，砂土與黏土的土-水特征曲線如圖2所示。

表1 二種非飽和土體的計算參數(shù)[20]

圖2 砂土與黏土的土-水特征曲線

3.2 物理力學參數(shù)

砂土與黏土的力學參數(shù)根據(jù)室內試驗確定，采用Mohr-Coulomb彈塑性準則，其物理力學性質參數(shù)[20]見表2。

表2 砂土與黏土的物理力學參數(shù)[20]

3.3 滲流各向異性定義及計算工況

以往理論、試驗及數(shù)值研究中大多忽略邊坡土體的各向異性程度及各向異性方向，然而事實上無論是黏土，砂土等邊坡土質均存在各向異性，對于滲透矩陣[C]而言，一般表達式可以表達為：

(4)

式(4)中C11=kxcos2α+kysin2α，C22=kxsin2α+kycos2α，C12=C21=kxsinαcosα+kysinαcosα。其中，kx、ky、α根據(jù)圖2定義，kx為水平向滲透系數(shù)，ky為垂直向滲透系數(shù)，α為ky與y軸的方向夾角。當α=0°時，[C]退化為：

(5)

式(5)是各向異性滲透矩陣，定義各向異性程度kr=kx/ky，但是各向異性不僅僅存在水平豎直滲透系數(shù)不一致的問題，同時存在各向異性方向α的問題。

為綜合討論邊坡不同土質的各向異性，設置以下工況：滑坡體土質采用黏土及砂土，各向異性程度kr=1、10、50、100,同時各向異性方位角α=0°、15°、30°、45°、60°、75°及90°。庫水位從175 m高程變化到145 m高程，下降速率為1.2 m/d，同時考慮到庫水位下降停止后的55 d，總共計算時間為80 d。計算工況如表3所示。

表3 計算工況

4 計算結果與分析

根據(jù)表3的計算工況對砂土與黏土進行56組工況的計算，依照砂土及黏土分類分析各向異性程度kr變化對邊坡滲流穩(wěn)定的影響。

4.1 各向異性對滲流特性的影響

4.1.1 砂土

砂土邊坡不同各向異性程度kr及方向α下的浸潤線變化如圖3所示。由圖3可見：庫水位下降停止時刻，在接近庫岸部位的浸潤線存在上凸現(xiàn)象，即岸坡內部的浸潤線變化要滯后于庫區(qū)水位的變化。不同各向異性方向α對浸潤線變化影響巨大，隨著各向異性角度α的增大，庫水位下降停止時刻(第25 d)岸坡內部的浸潤線高程逐漸抬高，但是增加高程卻逐漸衰減，對于α大于45°情況下浸潤線幾乎一致。

圖3 砂土各向異性程度及方向對浸潤線變化影響

不同各向異性程度kr對浸潤線的影響也較大。隨著各向異性程度kr的增大，浸潤線抬高迅速，對于kr大于50的情況，浸潤線高程幾乎一致。

4.1.2 黏土

黏土邊坡不同各向異性程度kr及方向α下的浸潤線變化如圖4所示。由圖4可見：黏土邊坡的整體浸潤線比砂土邊坡更高，同時，不同各向異性程度及方向下的浸潤線差異較之砂土邊坡更小，各向異性程度kr越大，各向異性方向角α越大，浸潤線的高程越高。

圖4 黏土各向異性程度及方向對浸潤線變化影響

4.1.3 浸潤線滯后高程分析

通過4.1.1及4.1.2節(jié)可以發(fā)現(xiàn)，無論是砂土還是黏土邊坡，庫水位驟降下岸坡內部浸潤線變化均滯后于庫岸的水位，這種滯后性會導致局部的滲透坡降陡增，滲流力指向坡外，加劇岸坡失穩(wěn)。因此為定量化研究這種滯后效應，定義庫水位驟降下的浸潤線滯后高程為圖5a中下部監(jiān)測面與庫水位驟降停止(第25 d)的浸潤線的交點高程至死水位145 m高程的高度。圖5b、圖5c分別為不同滲流各向異性程度及方向下的砂土邊坡與黏土邊坡的浸潤線滯后高程。

總體上看，隨著各向異性程度kr及方向角α的增大，砂土邊坡與黏土邊坡的浸潤線滯后高程逐漸增大至極大值。這是因為各向異性方向角α增大，對于kr>1的情況而言，與x軸平行方向的滲透系數(shù)越來越小，坡體內部的水很難排滲到庫區(qū)，因此水位滯留現(xiàn)象越明顯，因此浸潤線滯后高程越大。同時，各向異性程度kr越大，土體的豎直向的滲透系數(shù)也越小，庫水位較難下滲，因而浸潤線滯后高程也越大。對比砂土與黏土邊坡，可以發(fā)現(xiàn)滲流各向異性特性對砂土邊坡的浸潤線滯后高程影響巨大，僅考慮各向異性程度kr而忽略各向異性方向α情況下最大與最小浸潤線滯后高程相差6.51%，而同時考慮各向異性程度kr及方向α情況下最大與最小浸潤線滯后高程相差35.53%。對于黏土邊坡，滲流各向異性特性影響較小，僅考慮各向異性程度kr而忽略各向異性方向α情況下最大與最小浸潤線滯后高程相差4.7%，而同時考慮各向異性程度kr及方向α情況下最大與最小浸潤線滯后高程相差9.94%。

圖5 不同各向異性程度及方向砂土與 黏土邊坡的浸潤線滯后高程

4.2 各向異性對滑坡不同位置變形的影響

4.2.1 砂土

基于孔壓-應力耦合方法，得到了砂土不同各向異性下滑坡不同監(jiān)測面(上部監(jiān)測面、中部監(jiān)測面及下部監(jiān)測面)的水平位移分布規(guī)律，如圖6所示。

由圖6可見:庫水位驟降停止時刻不同滑坡位置的位移變化各不相同，對于上部監(jiān)測面而言，水平向位移隨高程呈現(xiàn)先增大后減小的變化規(guī)律，其最大位移發(fā)生在距地表約為5 m處。而中部與下部監(jiān)測面的水平位移則是隨著高程逐漸減小。距離庫岸越近，監(jiān)測面的整體位移越大，結合4.1節(jié)分析可知，庫水位下降靠近庫岸岸坡內部水位線滯后，對周圍土體骨架產生向臨空面及向下的滲流動水壓力(動態(tài)擴張力)，牽引帶動坡面土體[6-8]，因此可以判定庫水位下降坡趾最容易發(fā)生破壞。不同各向異性特性對庫水位停止時的水平位移影響較大。上部監(jiān)測面位移數(shù)值差異較小，各向異性方向角度α越大，下部監(jiān)測面的整體位移越大，但是中部監(jiān)測面的位移卻越?。桓飨虍愋猿潭萲r越大，中部監(jiān)測面水平位移越小，而下部監(jiān)測面的水平位移越大。這可能是因為各向異性方位角α越大，各向異性程度kr越大，坡趾位置的滲透坡降逐漸增大，而坡中位置的滲透坡降逐漸減小所致。

圖6 不同各向異性特性下砂土水平位移變化

4.2.2 黏土

黏土邊坡不同各向異性特性對水平位移影響規(guī)律如圖7所示。

圖7 不同各向異性特性下黏土水平位移變化

由圖7可見：不同各向異性特性對黏土邊坡的水平位移影響規(guī)律與砂土邊坡類似，但是值得注意的是，各向異性程度kr及各向異性方向α對黏土邊坡的影響比砂土邊坡影響要小，體現(xiàn)在不同各向異性程度kr及各向異性方向α下的水平位移差異較小。與砂土邊坡對比可知，黏土邊坡的上部監(jiān)測面與中部監(jiān)測面的水平位移值整體上均小于砂土邊坡，而下部監(jiān)測面的水平位移值卻大于砂土邊坡，這也是容易理解的，因為黏土邊坡的坡降突變集中于坡趾，庫水位下降對坡中及坡頂?shù)挠绊懖淮螅鴰焖惑E降不僅僅影響到砂土邊坡的坡趾，對坡中還具有一定的影響。

4.2.3 最大水平位移分析

為定量化研究不同各向異性特性下砂土與黏土邊坡不同位移的最大水平位移，對不同工況的最大水平位移進行了統(tǒng)計，其統(tǒng)計結果如圖8所示。

圖8 不同各向異性特性下砂土與黏土的最大水平位移。

由圖8可見：不同各向異性特性對最大水平位移的影響隨監(jiān)測面位置的不同而不同。下部監(jiān)測面的位移差異最大，中部監(jiān)測面次之，而上部監(jiān)測面則最小。僅考慮各向異性程度kr但忽略各向異性方向α下砂土上部、中部及下部監(jiān)測面的差異分別為69.5%、20.7%、6.9%，黏土上部、中部及下部監(jiān)測面的差異分別為75.05%、77.2%、21.4%；同時考慮各向異性程度kr及方向α情況下砂土上部、中部及下部監(jiān)測面的差異在80.0%、94.6%、45.0%，黏土則為75.1%、87.1%、32.0%。

4.2.4 庫水位下降對邊坡的滲流穩(wěn)定性的影響

圖9是各向異性程度kr=1，各向異性方向α=0°情況下砂土邊坡的庫水位動態(tài)變化圖。圖10顯示：在庫水位下降過程中，土體的持水性使浸潤線的下降滯后于庫水位，因此對邊坡的穩(wěn)定性影響主要體現(xiàn)在以下三個方面，一是坡體內具有較高的孔隙水壓(超孔隙水壓力)且不能及時消散，形成非穩(wěn)態(tài)滲流場，對周圍土體骨架產生向臨空面及向下的滲流動水壓力(動態(tài)擴張力)，同時由水頭差產生的靜水壓力相當于對前緣施加了一個推力，運動著的地下水“拖拽”土體或土粒向滲流方向前進，這是對邊坡穩(wěn)定的不利方面；二是庫水位在較高水平情況下庫區(qū)水壓力較高，對庫岸滑坡起到阻滑的作用，但是當庫水位下降至低水位時，上游水壓力卸載，使得阻滑力降低，此為邊坡穩(wěn)定的不利影響；三是庫水位降低導致岸坡內部的孔隙水壓力降低，土體的有效應力及強度提高，此為邊坡穩(wěn)定的有利影響。綜上所述，庫區(qū)岸坡穩(wěn)定性受以上三個因素的綜合影響，其中，邊坡土質、滲流各向異性特性則是影響邊坡穩(wěn)定的重要因素。

圖9 庫水位驟降不同時刻浸潤線變化

由本文數(shù)值模擬可知庫水位驟降下坡趾的水平位移最大，這同樣也是由于水位線“滯后”導致土體內部指向坡面外的拖拽力導致。本文基于孔壓-應力耦合方法得到庫水位驟降全過程的水平位移云圖(圖10a)顯示：滑坡坡趾處首先產生位移滑動，隨著庫水位的進一步下降，坡中及坡頂位移逐漸增大，可見庫水位下降滑坡破壞首先發(fā)生在坡趾。江強強等[21]對庫水位下降下的滑坡模型的室內試驗發(fā)現(xiàn)，庫水位驟降下滑坡的坡趾處首先產生貫穿裂縫，隨后發(fā)生局部流滑動最終發(fā)生滑坡前緣的整體破壞，這與本文的數(shù)值模擬的結論是一致的，從而也驗證了本文數(shù)值模擬及孔壓-應力耦合方法的正確與合理性。因此，在實際庫區(qū)滑坡治理工程中，應該加強滑坡坡趾處的加固措施。

圖10 本文數(shù)值模擬與文獻[21]試驗結果對比

4.3 各向異性對滑坡安全系數(shù)的影響

4.3.1 砂土

滲流各向異性對砂土邊坡的安全系數(shù)影響規(guī)律如圖11所示，由圖11可見：安全系數(shù)隨庫水位驟降呈現(xiàn)先減小后增大的規(guī)律(除α=45°情況下的kr=10及kr=100)，值得注意的一點是滲流各向異性特性不僅僅對邊坡的安全系數(shù)大小有影響，同時對最小安全系數(shù)的時刻具有影響。滲流各向異性角度α越大，各向異性程度kr越大，邊坡的整體安全系數(shù)越小，這與4.1節(jié)討論的庫水位下降下邊坡內部的浸潤線滯后效應有關。同時，滲流各向異性角度α及各向異性程度kr的增大同時也使得邊坡最小安全系數(shù)出現(xiàn)的時刻提前，這是由于庫水位下降情況下存在三個作用過程，一個是庫岸坡面的水壓力卸載，使得岸坡阻滑力降低，此為邊坡穩(wěn)定的不利因素，一個是岸坡內部浸潤線滯后，導致水力坡降增大，此亦為不利因素，而另一個則是庫岸內部的水位線降低，土體的有效應力及土體強度升高，此為邊坡穩(wěn)定的有利因素，結合4.1.1節(jié)分析結果可知，當滲流各向異性角度α及各向異性程度kr較小時，邊坡的排滲能力較強，此時庫岸坡面的水壓力卸載及水位線滯后的影響要小于庫岸內部的水位線降低導致土體的有效應力及土體強度升高的影響，隨著各向異性角度α及各向異性程度kr逐漸增大，邊坡的排滲能力降低，所以最小安全系數(shù)出現(xiàn)的時間逐漸滯后。

圖11 砂土滲流各向異性對安全系數(shù)的影響規(guī)律

4.3.2 黏土

滲流各向異性對黏土邊坡的安全系數(shù)影響規(guī)律如圖12所示，由圖12可見：各向異性程度及方向對黏土邊坡的影響規(guī)律與砂土邊坡類似，但是與砂土邊坡相比，滲流各向異性特性對黏土邊坡的影響更小。值得注意的是，黏土邊坡的安全系數(shù)整體上要比砂土邊坡更低，這也是由于黏土邊坡庫水位下降時的浸潤線滯后效應更加明顯導致。

圖12 黏土滲流各向異性對安全系數(shù)的影響規(guī)律

4.3.3 最小安全系數(shù)

滲流各向異性對砂土與黏土邊坡的最小安全系數(shù)影響規(guī)律(圖13)顯示：隨著各向異性程度kr及各向異性方向α的增大，最小安全系數(shù)逐漸走向最低。值得注意的是，黏土邊坡在庫水位驟降下的整體最小安全系數(shù)要小于砂土邊坡。僅考慮各向異性程度kr但忽略各向異性方向α下砂土邊坡的最小安全系數(shù)相差2.72%，黏土邊坡則相差9.29%，而綜合考慮各向異性程度kr及各向異性方向α砂土邊坡的最小安全系數(shù)相差11.58%，而黏土邊坡則為12.31%。

圖13 砂土與黏土邊坡的最小安全系數(shù)

5 結論

本文基于孔壓-應力耦合模型，以三峽庫區(qū)蔡坡堆積體為例，同時考慮到砂土質與黏土質邊坡的各向異性程度kr及各向異性方向α的影響，對庫水位下降下邊坡的滲流及穩(wěn)定性進行了數(shù)值模擬，得到了以下結論：

(1)各向異性程度kr及其方向α的增大降低了滑坡土體的排滲能力，在庫水位下降過程中岸坡內部的浸潤線出現(xiàn)“滯后”效應，使得邊坡內部孔壓較高，位移增大及安全系數(shù)降低。

(2)各向異性程度kr及方位角α越大，浸潤線滯后高程、最大水平位移增大，最小安全系數(shù)的值也越大，砂土邊坡的浸潤線滯后程度要小于黏土邊坡，但是最大水平位移和最小安全系數(shù)要大于黏土邊坡。

(3)對于砂土邊坡而言，僅考慮各向異性程度kr與同時考慮各向異性程度kr及方向α下的差異較大，而對于黏土邊坡兩者差異較小。因此對于滲透系數(shù)較大的砂土邊坡而言，必須全面考慮各向異性程度kr及其方向α的影響，但是對于滲透系數(shù)較小的黏土邊坡而言，可以適當忽略各向異性方向α以簡化計算。