降雨入滲下非飽和土質邊坡穩定性的最小勢能法研究

黃 強，程勇剛，周 偉，常曉林

(1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072； 2.武漢大學 水工巖石力學教育部重點實驗室，湖北 武漢 430072)

0 引 言

降雨是邊坡失穩的主要誘因之一[1,2]，其實質在于降雨入滲導致巖土體強度特性發生變化以及應力狀態調整，因此，研究一種可以考慮非穩定非飽和滲流作用的邊坡穩定性分析方法具有重要意義和工程實用價值。關于非穩定非飽和滲流下的邊坡穩定性分析方法國內外已經有許多研究。Fredlund[3]、黃潤秋等[4]分析了邊坡土體中基質吸力對穩定性的影響，認為降雨入滲導致基質吸力降低，土體強度減弱，從而使邊坡穩定性下降。Fredlund[5]比較了各種條分法的不同點，分析其在非飽和土質邊坡穩定性分析的問題。Alonso 等[6]考慮了土體類型、土體滲透性、土-水特征曲線、降雨持時和降雨強度等影響因素，利用極限平衡法進行二維非飽和土質邊坡的穩定性分析。Cai 等[7]、Huang 等[8]和榮冠等[9]將非穩定非飽和滲流有限元分析與有限元強度折減法相結合，對降雨或水位變化下的土坡穩定進行分析。李榮建等[10]利用極限平衡法和有限元強度折減法，分別探討了Fredlund 和Bishop 非飽和抗剪強度理論在邊坡穩定分析問題中的特點。

Mauldon and Ureta[11]從能量的角度出發，基于最小勢能原理，提出了一種針對棱柱形滑面的邊坡穩定性分析方法。李小強等[12]在此基礎上，利用最小勢能法推導了二維均質和非均質邊坡安全系數的計算公式，并與簡化Bishop 法結果進行對比。李鈾等[13]考慮滑床剪切勢能，進行錨固邊坡的穩定性分析。沈愛超等[14]將最小勢能法應用到任意滑移面的邊坡穩定分析中。溫樹杰等[15]研究了任意滑裂面的三維均質邊坡穩定分析的最小勢能法。

目前廣泛應用的條分法包括Mogenstern-Price法、Spencer法和Janbu法等，均是在土體摩爾-庫倫破壞準則以及力與力矩平衡條件理論下建立的。它們需要通過對土條側向力作不同假定使得問題靜定可解，同時計算過程中存在迭代收斂問題[16,17]。最小勢能法是基于整個滑體能量的理論體系，適用于任意形狀的滑面，不需要劃分土條，也就不存在條間力的假定以及迭代計算，同時根據最小勢能原理滿足平衡微分方程和力的邊界條件，具有廣闊的應用前景。現有的研究工作表明了最小勢能法在邊坡穩定分析問題中的有效性和可靠性，但沒有考慮非飽和滲流對邊坡穩定的作用，而基質吸力對土體抗剪強度有著重要影響[18,19]。

為此，本文引入不同時刻滑面上的孔隙水壓力來考慮非穩定非飽和滲流對邊坡系統勢能的作用，建立勢能函數，并基于Fredlund 雙應力變量理論將基質吸力的時空分布轉換為邊坡土體強度參數的時空分布，從而推導出非穩定非飽和滲流條件下的邊坡安全系數計算公式。最后通過算例，與現有方法進行對比，驗證了本文方法的合理性。

1 降雨入滲下非飽和土質邊坡穩定性的最小勢能法分析理論

1.1 均質土坡穩定性分析的最小勢能法

對于如圖1所示的單一均質邊坡ACBD，重度為γ，黏聚力為c′，內摩擦角為φ′。對于任意形狀滑面ADB，曲線方程可表示為y=f(x)的函數形式，點D為浸潤線與滑面交點。采用如下假定：

(1)滑體ACBD為剛性體，在滑裂面上發生彈性變形；

(2)對于滑裂面上任一微分段dl，其彈性變形按法向彈簧模擬，其剛度k與長度dl成正比；

(3)滑體ACBD整體在合外力R作用下發生一微小位移d，使得系統總勢能最小。

圖1 單一均質邊坡穩定分析示意圖Fig.1 stability analysis of single uniform soil slope

根據上述條件可也確定最小勢能法邊坡穩定性分析中的一些基本幾何與力學參數：

(1)

k=mdl

(3)

式中：dl為微分段長度；k為彈簧剛度；m為與材料有關的常數；n為外法線向量；n1，n2為在坐標軸的投影分量。

假設滑體為剛形體，在合外力作用下發生一微小位移，當系統恢復平衡狀態時，系統勢能最小，此時滑體在滑裂面上存儲的彈性勢能可表示為：

(4)

1.2 考慮非飽和非穩定滲流的最小勢能法

在涉及滲流的邊坡穩定性分析問題中，對于滑體的受力狀態需要考慮水對土體的作用力。以浸水土體為研究對象時，水對土顆粒的作用力可以看為兩個部分：①滑體容重的改變；②滑面上水的作用。對于第一部分，在計算土體自重時需要按飽和非飽和情況考慮。對于第二部分，在非飽和滲流問題中，存在滑面與自由面相交的情況。此時，位于水位線上下的滑面受孔隙水壓力作用，對系統做功；位于水位線以上的滑面，其處于非飽和區，土體抗剪強度受飽和度影響。

為此，在計算系統總勢能時需要考慮滑裂面上孔隙水壓力所做的功。

對于一些簡單的滲流問題，浸潤線可以用形如y=H(x)的函數表示，此時滑面上的孔隙水壓力分布為：

p=γ[H(x)-f(x)]

(5)

則滑面ADB段孔隙水壓力所做的功WH為：

(6)

式中：p為孔隙水壓力值；d1，d2分別為位移d在坐標軸的分量。

對于復雜的非飽和非穩定滲流問題，浸潤線很難以函數形式顯式表示，但當滑裂面上第i點的孔隙水壓力pi已知時，可對滑裂面上孔隙水壓力所做的功進行分段求和，WH近似為

(7)

則系統總勢能V為：

(8)

式中：Wτ為除滑裂面上孔隙水壓力外其他外力所做的功。

根據最小勢能原理，當系統勢能最小時，滿足：

(8)

求解式(9)即可得到位移分量d1，d2。對于單一均質土體，求解得到的d1，d2是以m為分母的表達式，即：

(10)

式中：U1(x)，U2(x)為不包含m項的函數表達式。

則滑面上的法向應力為：

f(σ)=m(d·n)=m[d1n1(x)+d2x2(x)]

(11)

根據法向應力分布，結合抗剪強度公式，沿滑面路徑進行積分即可確定滑體的總抗滑力。可以看出，最小勢能法是對整個滑體進行分析，在確定法向應力分布的過程中無需劃分條塊迭代，相比極限平衡法，計算過程簡單。

對于部分位于水位線以上的滑面，其處于非飽和區，滑面上土體的抗剪強度受基質吸力影響。對于這種非飽和土特性，采用Fredlund雙應力變量理論[20]，將基質吸力作為一個獨立的應力變量，通過納入到等效黏聚力中來考慮其對強度的影響。非飽和土抗剪強度公式為：

τf=c′+(σ-ua)tanφ′+(ua-uw)tanφb

(12)

式中：c′，φ′為有效黏聚力和有效內摩擦角；(σ-ua)為凈法向應力；(ua-uw)為基質吸力；τf為非飽和土抗剪強度；σ為總法向應力；ua為孔隙氣壓力；uw為孔隙水壓力；tanφb為與基質吸力有關的吸力內摩擦角。

Vanapalli and Fredlund[21]建立了φb與體積含水率之間的關系式：

(13)

式中：θw為體積含水率；θr為殘余體積含水率；θs為飽和體積含水率。

1.3 非穩定滲流條件下的最小勢能法邊坡安全系數定義

不考慮氣相的影響，將最小勢能法確定的法向應力分布代入非飽和抗剪強度公式中，從而獲得滑面上土體抗剪強度分布：

滑面上每一微分段在滑動方向上能提供的極限抗滑力為：

dT=cosα·τfdl

(16)

式中：dT為每一微分段能提供的極限抗滑力；α為抗滑力與滑動方向的夾角；t為外切線向量。

則總極限抗滑力為：

(17)

當滑體只受重力W作用時，在飽和非飽和滲流條件下，土體自重按飽和非飽和情況計算。此時下滑力為：

(18)

則最小勢能法安全系數為：

(19)

對于圓弧形滑面，總極限抗滑力的表達式為：

(20)

則安全系數按抗滑力矩與下滑力矩的比值定義：

(21)

式中：r為圓弧滑面半徑；rw為下滑力力臂。

2 降雨入滲下非飽和土質邊坡穩定性的最小勢能法實現過程

上述分析是針對某一指定潛在滑面的，在邊坡穩定性分析問題中，為確定最危險滑面的位置以及最小安全系數，需要選取若干個潛在滑面進行對比分析。因此，降雨入滲下非飽和土質邊坡的最小勢能法分析流程如圖2所示，其分析過程如下：

圖2 降雨入滲下非飽和土質邊坡的最小勢能法流程圖Fig.2 Flowchart of Minimum potential energy method for stability analysis of unsaturated soil slope under rainfall infiltration

(1)邊坡模型的建立以及土體參數確定；

(2)降雨入滲下的非飽和土質邊坡滲流分析，確定孔隙水壓力分布；

(3)擬定若干個可能的滑移面，獲取滑面上的孔隙水壓力及浸潤線位置，確定滑體在飽和非飽和情況下的自重，飽和區滑面上的水壓力及非飽和區的基質吸力；

(4)基于最小勢能原理，求解各滑體在滿足勢能最小原則時的位移，以確定滑面上的法向應力分布；

(5)分別計算各個滑體的抗滑力及下滑力，按式(19)或式(21)求解各個滑面的安全系數；

(6)對比選取安全系數最小的滑動面為最危險滑面，其安全系數為邊坡的最小安全系數。

3 基于最小勢能法的降雨入滲下非飽和土質邊坡穩定分析

為了簡化問題，下面分析沒有考慮滑面上的剪切勢能，但其影響不大[13]。

3.1 降雨入滲下的非飽和均質粉砂土邊坡穩定分析

取一簡單的降雨入滲模型[7]，按圓弧滑面考慮。堤高10 m，堤頂寬24 m，堤坡比1∶1.5，跛腳距模型左邊界16 m，模型總高度15 m，總寬度55 m。模型兩側為法向約束，底部為全約束，初始地下水位高度為5 m，模型BC與BD段受降雨作用，降雨強度為5 mm/h，降雨持續10 d。模型共劃分9 464個單元，9 706 個節點。計算模型如圖3所示，土體參數見表1，非飽和土-水特征曲線(見圖4)采用Van-Genuchten模型[22]。

3.1.1 降雨入滲下非飽和土質邊坡的滲流場分析

圖5給出了降雨不同時刻下邊坡土體的含水率分布情況。

圖3 計算模型及網格Fig.3 Calculation model and finite element mesh

圖4 土-水特征曲線Fig.4 Soil-water characteristic curve

表1 材料參數Tab.1 parameters of material

圖5 降雨入滲下邊坡土體不同時刻的體積含水率云圖Fig.5 Contours of water content at different times for a soil slope under rainfall conditions

在初始地下水位作用下，邊坡土體的體積含水率呈層狀分布。降雨開始后，坡頂和坡面迅速飽和。伴隨降雨繼續，表面雨水從暫態飽和區向深部移動，形成由坡頂往深部體積含水率飽和(高)→非飽和(低)→飽和(高)的封閉現象，并隨著降雨入滲深度增大，封閉區逐漸縮小且向內部發展。

圖6給出了坡頂C處斷面在降雨入滲條件下壓力水頭隨高程的變化情況。從圖6中可以看出，在初始水位下，壓力水頭隨高程呈線性分布，當降雨開始，坡頂表面迅速飽和。隨著降雨的繼續，濕潤鋒向下推進，表面飽和區不斷擴大，土體基質勢梯度逐漸減小，同時底部部分土體吸水飽和，并逐漸向上部發展，使得坡體中部的非飽和帶縮小。

圖6 不同時刻下壓力水頭沿高度的分布情況Fig.6 Development of water pressure head at different times for various heights of slope

3.1.2 不同方法的邊坡安全系數對比分析

分別按最小勢能法、Spencer法和Morgenstern-Price法對降雨入滲下的非飽和土質邊坡進行穩定性分析。圖7給出了3種方法確定的邊坡安全系數隨降雨歷時的演化過程，圖8給出了3種方法在不同降雨歷時下確定的最危險滑面位置，黑色虛線為Spencer法確定的最危險滑面，灰色虛線為Morgenstern-Price法確定的最危險滑面，黑色粗線為最小勢能法(MPE)確定的最危險滑面，點劃線為該時刻下的浸潤線。。

圖7 3種方法確定的安全系數演化過程Fig.7 factors of safety histories between minimum potential energy method, Spencer method and Morgenstern-Price method

圖8 3種方法確定的不同時刻下最危險滑面位置Fig.8 Most critical slip surfaces at different times between minimum potential energy method, Spencer method and Morgenstern-Price method

由圖7可知，最小勢能法、Spencer法和Morgenstern-Price法確定的安全系數變化規律基本相同，總體表現為安全系數隨降雨時間不斷減小。降雨開始時，雨水入滲，土壤中基質吸力減小，抗剪強度降低，安全系數隨之下降。在降雨大約持續3天后，雨水滲入坡體一定深度，安全系數迅速下降。隨著降雨的繼續，土體孔隙被水充填，基質吸力不斷減小，雨水入滲速率下降，使得安全系數下降趨勢逐漸放緩。同時比較3種方法在不同時刻下的邊坡安全系數，其基本一致，最小勢能法與Spencer法、Morgenstern-Price法最大相差0.054，說明本文提出的降雨入滲下非飽和土邊坡穩定性的最小勢能法是合理的。從圖8可以發現，3種方法確定的最危險滑面位置基本一致，Spencer法和Morgenstern-Price法確定的最危險位置完全一致，最小勢能法除降雨1天時的滑面位置相比Spencer法要淺一些，其他時刻的滑面位置相同，證明了最小勢能法在降雨入滲下非飽和土邊坡穩定性分析中的合理性。為了簡化問題，本文沒有考慮滑面上的剪切勢能，但其影響不大[13]。

3.2 降雨入滲下的非飽和均質黏土邊坡穩定分析

取一坡高5 m、坡長30 m、坡度為1∶2的典型均質黏土邊坡[23]。初始地下水位高度為3 m，模型左邊界為透水邊界，模型右邊界及底部不透水，模型ABCD整段受降雨作用，降雨強度為3.125 mm/h，降雨持續48 h。模型共劃分3 613個單元， 3 754個節點。計算模型如圖9所示，土體參數見表2，非飽和土-水特征曲線(見圖10)采用Van-Genuchten模型[24]。

圖9 計算模型及網格Fig.9 Calculation model and finite element mesh

表2 材料參數Tab.2 parameters of material

圖10 土-水特征曲線Fig.10 Soil-water characteristic curve

3.2.1 不同方法的邊坡安全系數對比分析

分別按最小勢能法、Spencer法和Morgenstern-Price法對降雨入滲下的非飽和土質邊坡進行穩定性分析。圖11給出了3種方法確定的邊坡安全系數隨降雨歷時的演化過程，圖12給出了3種方法在不同降雨歷時下確定的最危險滑面位置，黑色虛線為Spencer法確定的最危險滑面，灰色虛線為Morgenstern-Price法確定的最危險滑面，黑色粗線為最小勢能法(MPE)確定的最危險滑面，點劃線為該時刻下的浸潤線。

圖11 3種方法確定的安全系數演化過程Fig.11 factors of safety histories between minimum potential energy method, Spencer method and Morgenstern-Price method

圖12 3種方法確定的不同時刻下最危險滑面位置Fig.12 Most critical slip surfaces at different times between minimum potential energy method, Spencer method and Morgenstern-Price method

由圖11可知，最小勢能法、Spencer法和Morgenstern-Price法確定的安全系數變化規律基本相同，總體表現為安全系數隨降雨時間不斷減小，降雨結束后安全系數逐漸回升。降雨開始時，雨水入滲，土壤中基質吸力減小，抗剪強度降低，安全系數隨之下降。降雨停止后，孔隙水在自重作用下消散，基質吸力增大，抗剪強度逐漸增加，安全系數回升。同時比較3種方法在不同時刻下的邊坡安全系數，其基本一致，最小勢能法相比Spencer法、Morgenstern-Price法安全系數最大相差0.030，說明本文提出的降雨入滲下非飽和土邊坡穩定性的最小勢能法是合理的，同時也說明傳統非飽和土質邊坡穩定性分析方法偏保守。從圖12可以發現，3種方法確定的最危險滑面位置基本一致，最小勢能法的滑面位置較其他兩種方法要淺一些，但總體相差不大。這進一步證明了最小勢能法在降雨入滲下非飽和土邊坡穩定性分析中的合理性。

4 結 語

本文通過建立考慮滲流作用的勢能函數，基于Fredlund非飽和土抗剪強度理論，提出了非穩定非飽和滲流下邊坡穩定性的最小勢能法，并對降雨入滲下的邊坡穩定性進行分析，得出如下結論：

(1)考慮非飽和非穩定滲流的最小勢能法是將滑體視為一個整體，不假定滑面形狀，基于最小勢能原理，進行邊坡穩定性分析。該方法無需劃分條塊及迭代，不對條間力作假定，滿足平衡微分方程和力的邊界條件，可以方便納入慣性力及地面荷載等外部荷載，同時計算過程簡單，具有較好的實用價值。其適用于任意滑面的飽和-非飽和單一均質邊坡穩定分析問題，對于多層土的應用有待進一步研究。

(2)以降雨入滲下的邊坡穩定分析為例，對比Spencer法、Morgenstern-Price法和最小勢能法結果發現，三者確定的最危險滑面位置基本一致，且安全系數隨降雨時間的演化規律也基本相同，同時最小勢能法與其他兩種方法的安全系數最大相差0.054，驗證了本文方法的合理性，同時說明了傳統非飽和土質邊坡穩定性分析方法偏于保守。本文分析沒有考慮滑面上的剪切勢能，但其對結果影響不大。

□

[1] 黃潤秋. 20世紀以來中國的大型滑坡及其發生機制[J]. 巖石力學與工程學報, 2007,26(3):433-454.

[2] 周創兵, 李典慶. 暴雨誘發滑坡致災機理與減災方法研究進展[J]. 地球科學進展, 2009,24(5):477-487.

[3] Fredlund D G, Rahardjo H. Soil mechanics for unsaturated soils[M]. New York: Wiley, 1993.

[4] 黃潤秋, 戚國慶. 非飽和滲流基質吸力對邊坡穩定性的影響[J]. 工程地質學報, 2002,10(4):343-348.

[5] Fredlund D G, Krahn J. Comparison of slope stability methods of analysis[J]. Canadian Geotechnical Journal, 1977,14(3):429-439.

[6] Alonso E, Gens A, Lioret A, et al. Effect of rain infiltration on the stability of slopes[J]. Unsaturated Soils, 1995,1:241-249.

[7] Cai F, Ugai K. Numerical analysis of rainfall effects on slope stability[J]. International Journal of Geomechanics, ASCE, 2004,4(2):69-78.

[8] HUANG Mao-song, JIA Cang-qin. Strength reduction FEM in stability analysis of soil slopes subjected to transient unsaturated seepage[J]. Computers and Geotechnics, 2009,36(1-2):93-101.

[9] 榮 冠, 王思敬, 王恩志,等. 強降雨下元磨公路典型工程邊坡穩定性研究[J]. 巖石力學與工程學報, 2008, 27(4):704-711.

[10] 李榮建, 于玉貞, 鄧麗軍, 等. 非飽和土邊坡穩定分析方法探討[J]. 巖土力學, 2007,28(10):2 060-2 064.

[11] Mauldon M，Ureta J. Stability of rock wedges with multiple sliding surface[J]. Candian Geotechnital Journal, 1996,14(1):51-66.

[12] 李小強, 白世偉, 李 鈾. 最小勢能方法在二維邊坡穩定分析中的應用[J]. 巖土力學, 2004,25(6):909-912.

[13] 李 鈾, 陸 洋, 李 鈮. 錨桿(索)加固邊坡的最小勢能穩定分析方法研究[J]. 巖土力學, 2008,29(9):2 329-2 334.

[14] 沈愛超, 李 鈾. 單一地層任意滑移面的最小勢能邊坡穩定性分析方法[J]. 巖土力學, 2009,30(8):2 463-2 466.

[15] 溫樹杰, 羅 惠, 李 鈾, 等. 基于最小勢能原理的三維邊坡穩定性分析[J]. 巖石力學與工程學報, 2015,34(增1):3 298-3 305.

[16] Morgenstern N R, Price V E. The analysis of the stability of General Slip Surface[J]. Geotechnique, 1965,15(1):79-93.

[17] 陳祖煜. 土坡穩定分析通用條分法及其改進[J]. 巖土工程學報, 1983,5(4):11-27.

[18] 吳宏偉, 陳守義, 龐宇威,等. 雨水入滲對非飽和土坡穩定性影響的參數研究[J]. 巖土力學, 1999,20(1):1-14.

[19] 吳俊杰, 王成華, 李廣信. 非飽和土基質吸力對邊坡穩定的影響[J]. 巖土力學, 2004,25(5):732-736.

[20] Fredlund D G, Morgenstern N R, Widger R A. The shear strength of unsaturated soils[J]. Canadian Geotechnical Journal, 1978,15(3):313-321.

[21] Vannpalli S K, Fredlund D G, Pufahl D E, et al. Model for the prediction of shear strength with respect to soil suction[J]. Canadian Geotechnical Journal, 1996,33(3):379-392.

[22] VAN G M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(5):892-898.

[23] 蔣水華, 李典慶, 周創兵,等. 考慮參數空間變異性的非飽和土坡可靠度分析[J]. 巖土力學, 2014,35(9):2 569-2 578.

[24] 吳煜禾, 張洪江, 王 偉,等. 重慶四面山不同土地利用方式土壤水分特征曲線測定與評價[J]. 西南大學學報(自然科學版), 2011,33(5):102-108.

[25] Geo-Slope International Ltd. User's guide for finite element seepage analysis[M]. Canada: Geo-Slope International Ltd., 2007.

[26] Geo-Slope International Ltd. User's guide for slope stability analysis[M]. Canada: Geo-Slope International Ltd., 2007.