地震作用下基于嚴格極限平衡法的邊坡穩定性分析

吳敏強，蔡宜洲，周 豪

(三峽大學 水利與環境學院，湖北 宜昌 443002)

1 概 述

極限平衡法于1927年由Fellenius首次提出，是邊坡穩定性分析中最常用的方法[1]。當今應用較為廣泛的極限平衡法有Fellenius法[2]、簡化Bishop法[3]、Janbu法[4]、Sarma法[5]、Morgenstern-Price法[6]等。極限平衡法的研究領域分布廣泛，條間力的假設與滑動面的搜索是其主要應用方向，而且在三維擴展方面也成果顯著。

關于剛體極限平衡法的研究，國內學界大致有如下幾種觀點。朱大勇等[7]通過假設一組合式為零的條間剪力分布，使整個滑體所有平衡條件得到滿足，驗證了簡化Bishop法計算公式中沒有條塊間剪力，是因某種組合式為零，而不是條間剪力實際為零的設想。鄧東平等[8]找到了一種用隨機角來搜索隨機滑動面、簡化Janbu計算安全系數的新方法，通過隨機模擬曲線，發現隨機角物理意義明確，優勢獨特。鄭穎人等[9]針對Sarma法條間剪力方程存在的兩個問題：不滿足邊坡穩定性分析的合理性要求與不能正確表示任意條塊的條間剪力，給出了新的條間剪力方程。

在整個發展過程中，學者們一直將多余未知量的假設選在條間作用力或作用位置上。對于嚴格的二維極限平衡法而言，這種處理方法尚能夠得到令人比較滿意的安全系數結果，但有時也存在收斂性問題。對于三維極限平衡法而言，由于假設較多，例如假設合理性、力學嚴格性，這種研究思路已經漸顯不足。

Bell、朱大勇和鄭宏等學者將研究視線轉到滑面正應力的分布模式上來。這種極限平衡分析方法具有力學嚴格性，與其他嚴格的條分法相比，求解過程更便捷、易于編程實現且不存在收斂性問題。

嚴格極限平衡法發展潛力巨大，如采用此種方法對地震力作用下的邊坡穩定性進行研究，將會產生可靠有效的理論成果，也有望取得不錯的實際效益。如若研究成果的合理性與適用性被證明，并被廣泛接受，大量應用，便可在工程邊坡的穩定性分析中，獲得省時省力、節約成本、提升速率的良好效果。

2 基本原理

在滑面正應力修正法計算分析中，若能準確構建滑面上的正應力分布模式，就可以求解出準確的安全系數。對于實際計算過程，滑面正應力分布基本無法直接獲得。但若是假設一個合理的分布模式模擬真實狀態，計算精度達到要求也是允許的。對于滑面正應力的分布模式，研究工作主要存在于兩個部分：①滑面正應力假設模式；②求解安全系數。

2.1 滑面正應力假設模式

Bell提出將滑面正應力函數假設成含有兩個待定參數的方程，這樣就與安全系數這個未知數一起構成3個未知數，3個未知數3個方程理論上一定能求出數值解，這樣安全系數在該假設模式下一定可以求出。假設形式見式(1)。

(1)

式中：λ1，λ2為假設的未知參數；σ(x)為滑面正應力分布函數；σ0(x)為瑞典法的分布函數；a，b為滑體兩端點的橫坐標。

鄭宏教授提議將滑體看作積分域，通過數學原理將其轉變為邊界積分的研究方法不再局限于對滑體進行條分，可直接獲得其安全系數。他建立了條塊的平衡條件后發現，滑面正應力由兩部分構成，其中滑體體積力和條間作用力各占一部分。通過進一步的研究使滑面正應力的具體力學組成結構關系更加清晰明朗[10]。同時，在他的研究中提出了一種新的滑面正應力分布的假設模式：

(2)

2.2 求解安全系數

圖1為具有一般形狀滑動面的邊坡，y=s(x)為任意形狀滑動面；y=g(x)為坡表地形線；單位寬度條塊上重力w(x)；單位寬度條塊上水平地震力為kcw(x)；表面水平方向上的分布荷載為qx(x)、豎直方向上的分布荷載為qy(x)；滑面上孔隙水壓力為u(x)；總正應力為σ(x)；剪應力為τ(x)；y=yt(x)為邊坡體推力線；滑面形心點坐標為(xc,yc)。向坡體設置預應力錨索并視為一個外部集中力對邊坡整體進行極限平衡分析，其中P為每根錨桿(預應力錨索)所提供的錨固力，n為錨桿(預應力錨索)的數量，m為錨桿(預應力錨索)編號，i為錨桿(預應力錨索)在錨固過程中與水平線的夾角，(xpm,ypm)為錨桿形心點坐標。

圖1 邊坡滑體及受力分布

設邊坡安全系數為Fs，滑面摩擦系數為tanψ(x)，有效黏聚力為c(x)，根據摩爾-庫侖準則，可得：

(3)

安全系數的平衡方程如下：

(4)

(5)

(6)

式中：s′=ds/dx=tanα;α為計算點滑面的傾角。若能確定滑面正應力的分布函數，就可得到安全系數。

對邊坡滑體的任一條塊進行受力分析，見圖2。

圖2 邊坡滑體受力模型

圖2中，E為條間水平向作用力；T為條間豎向作用力；h為條塊高度；ht為條間力的作用點高度；Δx為條塊寬度。

(7)

(8)

根據摩爾-庫侖強度準則：

(9)

由式(9)可見，滑面正應力由滑體體積力和坡面外力對滑面正應力的貢獻(記為σ0(x))和土條間作用力對滑面正應力的貢獻(記為σ1(x))兩個部分構成，則可簡記為：

σ(x)=σ0(x)+σ1(x)

(10)

(11)

(12)

在滑面已知前提下，σ0(x)屬于已知函數。為了方程組能夠求解，可以采用以下形式來對滑面正應力進行假設：

σ(x)=σ0(x)+f(x)

(13)

式中：f(x)為逼近函數，f(x)構造為含有兩個待定參數的函數才能保證平衡方程能夠順利求解。逼近函數對滑面正應力的貢獻是非常有限的，不需過分追求與真實分布之間的吻合度，可以選用線性函數逼近：

(14)

將式(11)、式(13)、式(14)和摩爾-庫侖準則帶入平衡方程，經整理可得：

(15)

(16)

(17)

解平衡方程，可以得到：

(18)

其中，

(19)

(20)

(21)

其中，

T0=A3B2-A2B3

(22)

S0=A1B3-A3B1

(23)

G0=A1B2-A2B1

(24)

則式(18)可以寫成：

(25)

(26)

解式(25)，取實根得：

(27)

方程具有解析意義，是邊坡嚴格極限平衡法顯式解答，適用一般形狀滑裂面。

3 地震力作用下邊坡穩定性計算

3.1 均質邊坡穩定性計算

取某均質土質邊坡，邊坡長100 m，邊坡高20 m，邊坡坡角為45°，邊坡的斷面尺寸見圖3。該坡體重度為19 kN/m3，黏聚力為30 kPa，內摩擦角為35°。

圖3 土質邊坡斷面

均質邊坡的穩定性分別采用Slide軟件和Geostudio軟件來計算。由于邊坡的坡比為1∶1，可知該邊坡最大水平動力系數為2，最大豎向動力系數為1，地震烈度為Ⅶ度，參照《水利水電工程邊坡設計規范》(SL 386-2007)取地震作用綜合系數為 0.25，取地震作用重要性系數為 1.0。Slide軟件Bishop法計算得到安全系數為1.571，Janbu法計算得到安全系數為1.488，見圖4。

3.1.1 數值應力場獲取

通過對均質邊坡的二維建模以及材料參數的賦值，進行分析計算，以便從計算軟件中提取出滑動面上各節點的有用應力信息。土坡的穩定性計算分析采用理正巖土計算軟件，采用軟件中內置的地震烈度考慮地震效應。

圖4 均質邊坡穩定性計算

3.1.2 有限元法與滑面正應力修正法結合

初始正應力分布函數σ0(x)通過軟件得取可按照以下步驟操作來進行：

1) 條分整個滑體，讀取滑面穿過的條塊相對應位置的坐標。

2) 讀取相應坐標點的應力信息(σx,σy,τxy)，參照式(28)可以換算成對應點的初始正應力信息：

(28)

式中：條塊的底面中點與數值方向的夾角為β。單元的應力信息可直接由計算軟件導出,當正應力初始分布形式通過σ0(x)數值方法得到后，再結合線性修正方法，對任意邊坡安全系數的求解會更加方便、準確。

在Slide軟件中，根據斷面圖所提供的信息建立計算模型，由提供的已知條件賦予相應的荷載，地震烈度采用7級，可由計算軟件根據極限平衡理論自動搜索安全系數最小的滑面。在ANDIA軟件中提取滑面上的應力信息，經過公式轉換擬合出對應的初始滑面正應力分布函數，見圖5。

圖5 初始滑面正應力函數分布曲線

將得到的初始滑面正應力分布函數進行修正，計算得出邊坡在地震力作用下嚴格極限平衡法的安全系數。

經過曲線擬合，邊坡圓弧滑面正應力分布函數為：

y=0.0165x3-1.0566x2+17.38x+5.1878

計算在地震工況下邊坡的穩定性，將得到的滑面正應力的初始正應力分布函數σ0(x)，根據式(29)對滑面正應力進行修正，其中λ和ξ為假定未知參數。

σ(x)=σ0(x)+λζ1(x)+ξζ2(x)

(29)

(30)

(31)

修正后的滑面正應力分布函數見圖6。

圖6 修正后滑面正應力函數分布曲線

經過曲線擬合，修正后邊坡圓弧滑面正應力分布函數為：

y=0.0177x3-1.1266x2+15.366x+2.405

根據式(29)對滑面正應力進行修正，根據式(3)可以得出安全系數的顯示解。選取Geostudio軟件與嚴格極限平衡法結果進行比較，具體結果見表1。

表1 均質邊坡安全系數匯總

根據計算結果可以看出，基于極限平衡理論的滑面正應力修正法與數值方法相結合的穩定性計算，可以直接從數值模擬軟件中了解滑面底部的正應力分布情況，通過簡單的數值計算即可獲得比較準確的初始正應力，再對初始正應力進行修正，所得滑面正應力函數比較接近實際情況，所得安全系數的值與傳統極限平衡法所得結果的相對誤差非常小，從而驗證了嚴格極限平衡法的合理性，便于實際應用。

3.2 非均質邊坡穩定性計算

選取具有一般形狀典型的非均質邊坡，其剖面形狀見圖7。邊坡坡高15 m，邊坡比為1∶2，共有4個土層，各土層物理力學參數見表2。取地震峰值加速度為1.63 m/s2，對照中國地震烈度表取對應的地震烈度為Ⅶ度；參照《水利水電工程邊坡設計規范》(SL 386-2007)取地震作用綜合系數為 0.25，取地震作用重要性系數為1.0。

圖7 非均質邊坡剖面圖

表2 土層物理力學參數

基于數值應力場的非均質邊坡穩定性計算，可參照以下步驟進行：

1) 在Geostudio軟件中，根據剖面圖所提供的信息建立計算模型。

2) 可由計算軟件根據極限平衡理論自動搜索出安全系數最小的圓弧滑面。

3) 在大型有限元計算軟件ANSYS中，建立該剖面的二維有限元計算模型，其二維模型見圖8，采用D-P彈塑性本構模型。

4) 在有限元軟件中提取兩種滑動面上的應力信息，經過公式的轉換擬合出對應的初始滑面正應力分布函數。

5) 將得到的初始滑面正應力分布函數進行修正，計算得出安全系數。

6) 對比嚴格極限平衡法與傳統極限平衡法所得結果。

圖8 非均質邊坡二維有限元模型

在進行材料參數賦值后對該邊坡進行穩定性分析，對應極限平衡法搜索出的滑面位置在相應位置導出其應力信息，通過式(28)的轉換，可得到選取單元上的初始正應力值。經過正應力修正后便可以擬合出修正后滑面正應力分布函數，見圖9。

經過曲線擬合，修正后邊坡圓弧滑面正應力分布函數為：

y=-0.0071x3+0.2277x2+8.5385x+41.095

按照均質邊坡的計算過程，可以得出非均質邊坡在地震作用下嚴格極限平衡法的安全系數，其結果匯總見表3。

圖9 修正后滑面正應力函數分布曲線

表3 非均質邊坡安全系數匯總

根據計算結果可以看出，在地震作用下嚴格極限平衡法計算得出的結果小于Bishop法、M-P法和瑞典圓弧法的安全系數，但是大于Janbu法和Bishop法折線滑面的計算結果。與Bishop法相比，結果相對誤差為1.9%；與Bishop法折線滑面相比，結果相對誤差為2%；但是與Janbu法相比，相對誤差為11.1%?？傮w來說，誤差范圍在1.9%～11.1%之間。造成誤差的原因主要在于正應力提取的不完整和正應力修正造成一定的誤差，但在容許的范圍內，說明嚴格極限平衡法在非均質邊坡應用的合理性與適用性。

4 工程應用

4.1 工程概況

牙根水電站在雅礱江干流上，位于四川省甘孜州雅江縣境內，是我國重要水電能源基地雅礱江流域的大型電站工程?？碧浇衣秹沃穮^右岸Ⅰ-Ⅴ勘探線之間發育拉裂松動巖體，其地表后緣邊界高程自Ⅰ-Ⅴ線從2 640 m→2 667 m→2 690 m→2 697 m→2 597 m呈弧形變化，上游從PD17硐下部基巖突出處延伸至江邊，下游側從PD15硐延伸至江邊。此外，在Ⅶ線高高程PD31平硐(硐口高程2 612.67 m)硐深0～57 m，Ⅸ線低高程PD33平硐(硐口高程2 513.74 m)硐深0～129.5 m，中高程PD35平硐(硐口高程2 564.66m)等處亦發現有拉裂松動巖體分布。通過勘探，壩址區右岸存在的拉裂松動巖體對大壩的整體布局、溢洪道等水工建筑物全面布置，以及對施工和運行期的安全都存在隱患，評價其在不同工況下的穩定性，提出合適的加固處理措施。選取2號拉裂松動巖體為對象進行穩定性研究。

2號拉裂松動巖體分布于Ⅲ區，即Ⅶ-Ⅸ線下游側，拉裂松動巖體自Ⅶ線高高程斜向下延伸至Ⅸ線下游低高程江邊，呈現出從上游至下游自高向低的斜河向形態分布，橫向上受NE向陡傾角斷層的控制，其地表出露面積約6.3×104m2，體積約為300×104m3。該拉裂松動巖體內地表發育的斷層有：N20-40°E/SE∠80-85°的近橫斜河陡傾的f11-f14共4條Ⅲ級結構斷面(圖10)。

圖10 2號拉裂松動巖體平面地質

選取橫Ⅶ剖面對2號拉裂松動巖體的穩定性作計算。其計算剖面簡圖見圖11。

各地層物理力學參數見表4。

圖11 2號拉裂松動巖體橫Ⅶ剖面邊坡計算簡圖

表4 巖體物理力學參數表

4.2 計算成果分析

利用嚴格極限平衡法，計算邊坡在天然、暴雨、地震及地震+暴雨4種工況下的安全系數。暴雨工況表層參數c，f折減至90%，容重取飽和容重，地震工況的加速度為0.173 g，采用擬靜力法來計算，利用Slide軟件自帶的Bishop法和Janbu法計算不同工況下滑面的安全系數，結果見表5。

表5 計算結果匯總

注：表5中數值表示差值相對于Janbu法結果的百分比。

按照前述的方法計算不同工況下滑面的應力分布函數，按照式(29)進行修正并進行擬合，結果見圖12-圖15。

天然工況下圓弧滑面正應力分布見圖12。

則天然狀態圓弧滑面正應力分布函數為：

y=0.00008x3-0.051x2+8.1193x+187.36

帶入平衡方程中，解出該工況下的安全系數為1.193。

暴雨工況下圓弧滑面正應力分布見圖13。

則暴雨工況下圓弧滑面正應力分布函數為：

y=0.00007x3-0.0528x2+9.3091x+64.77

帶入平衡方程中，解出該工況下的安全系數為1.126。

地震工況下圓弧滑面正應力分布見圖14。

則地震工況下圓弧滑面正應力分布函數為：

y=0.0001x3-0.0678x2+10.537x+45.782

帶入平衡方程中，解出該工況下的安全系數為0.972。

地震+暴雨工況下圓弧滑面正應力分布見圖15。

則地震+暴雨工況下圓弧滑面正應力分布函數為：

y=0.00006x3-0.0481x2+8.9016x+53.2541

帶入平衡方程中，解出該工況下的安全系數為0.897。

圖12 天然工況下圓弧滑面正應力分布

圖13 暴雨工況下圓弧滑面正應力分布

圖14 地震工況下圓弧滑面正應力分布

圖15 地震+暴雨工況下圓弧滑面正應力分布

不同方法、不同工況計算結果匯總見表5。

通過對牙根水電站2號拉裂松動巖體采用基于數值應力場的滑面正應力修正法，分析邊坡在不同工況下的穩定性以及得出相對應的安全系數，并與Bishop法、Janbu法的安全系數作對比，可得出如下結論：3種極限平衡方法的計算結果相差在3%以內，三者得出的結論基本吻合。將滑面正應力修正法與有限元法相結合，理論明確，操作簡單，計算結果精度較高，適合大規模推廣。

5 結 語

本文采用滑面正應力修正法與數值應力場的分析方法對邊坡進行穩定性分析，其結果與傳統剛體極限平衡法的分析結果相比較相對誤差較小，進一步證明該方法可行性。通過對不同工況下的邊坡進行穩定性計算發現，邊坡在天然狀態下安全系數較小，當遇到暴雨或地震時安全系數下降幅度大，需要及時治理。通過對牙根水電站2號拉裂松動巖體采用基于數值應力場的滑面正應力修正法，分析邊坡在不同工況下的穩定性以及得出相對應的安全系數，并與Bishop法、Janbu法的安全系數作對比，可得出結論：3種極限平衡方法的計算結果相差在3%以內，三者得出的結論基本吻合。將滑面正應力修正法與有限元法相結合，理論明確，操作簡單，計算結果精度較高，適合大規模推廣。