邊坡穩(wěn)定分析的虛功率法1)

吳夢喜 楊家修 湛正剛

?(中國科學(xué)院力學(xué)研究所,北京 100190)

?(中國科學(xué)院大學(xué),北京 100049)

??(中國電建集團(tuán)貴陽水利水電勘測設(shè)計研究院,貴陽 550081)

引言

邊坡穩(wěn)定是巖土工程中的基本問題.由于大量自然和工程邊坡的穩(wěn)定安全裕度不大,近百年來巖土力學(xué)研究者與巖土工程師們持續(xù)不斷地努力發(fā)展邊坡穩(wěn)定分析方法,以提高邊坡穩(wěn)定分析的準(zhǔn)確性和便利性.邊坡穩(wěn)定分析方法可以分成剛體極限平衡法、有限元極限分析法、基于有限元應(yīng)力計算的改進(jìn)極限平衡法和極限分析法等.

極限平衡法通常也稱作條分法,將沿某一具有滑動可能的面上部的土體,切成若干有限寬度的豎條或斜條,把土條當(dāng)成剛體,根據(jù)靜力平衡條件求得滑動面上的滑動力和抗滑力,并以此為基礎(chǔ)確定邊坡的穩(wěn)定安全系數(shù).條分法是評估給定滑動面安全性的有力方法.最初出現(xiàn)的極限平衡法一般僅滿足土條的力和力矩平衡中的部分條件.如最早出現(xiàn)的瑞典條分法[1]僅滿足力矩平衡但不滿足力平衡、簡化Bishop 法[2]不滿足水平力的平衡、Janbu 簡化法[3]不滿足力矩平衡.這些方法的滑動面為圓弧,無需迭代就能求得安全系數(shù).因其簡便,目前仍被巖土工程師所采用.滑動面可以為任意形狀且同時滿足所有土條的力和力矩平衡的方法稱為嚴(yán)格的極限平衡法.嚴(yán)格的極限平衡法中,基于摩爾庫倫抗剪強(qiáng)度理論,對所有土體的強(qiáng)度除以同一個折減系數(shù),使每個土條底部的剪應(yīng)力與抗剪強(qiáng)度達(dá)到相等的臨界狀態(tài),折減系數(shù)即為給定滑動面的安全系數(shù).嚴(yán)格的極限平衡法求解的是靜不定問題[4],需要增加額外的假定.一般假定土條底部的法向力位于土條底部中點、且假定土條間法向力的作用位置或土條間作用力的方向后,變成方程比未知量多1 個的超靜定問題,安全系數(shù)需要迭代求解.基于不同的條間力的假定,嚴(yán)格的極限平衡法有十多個,如Spencer 法[5]、Morgenstern-Price 法[6]、一般極限平衡法(GLE)、Sarma 法[7],Janbu精確法.由于條間力的假定不同,不同的嚴(yán)格極限平衡法的計算結(jié)果存在差異,差異的大小與邊坡情況和所給定的滑動面位置有關(guān)[8].嚴(yán)格的極限平衡法,尚存在以下問題:(1)復(fù)雜的情況下難以找到比較符合實際的最危險滑動面; (2)求解的土條間的作用力為合力,滲流影響僅能通過浸潤線位置反映,土條間的作用力一般也不滿足摩爾庫倫強(qiáng)度準(zhǔn)則,因而,對于復(fù)雜滲流場和多層土邊坡情況,邊坡穩(wěn)定安全性評估結(jié)果往往與實際差異較大.

有限元極限分析法分為載荷增加法[8]和強(qiáng)度折減法[9-13],基于材料的彈塑性本構(gòu)關(guān)系,直接采用有限元求解邊坡處于塑性極限狀態(tài)時的載荷或材料強(qiáng)度折減系數(shù).此兩方法并無實質(zhì)差別.在邊坡穩(wěn)定分析中因強(qiáng)度折減法在已有的程序上實施更簡單而普遍采用.1975 年,Zienkiewicz 等[10]在研究土力學(xué)中的相關(guān)聯(lián)流動法則與非相關(guān)聯(lián)流動法則的文章中,在算例部分用有限元法分析了一個均勻邊坡的穩(wěn)定性,發(fā)現(xiàn)強(qiáng)度折減系數(shù)的倒數(shù)與極限平衡法計算的安全系數(shù)非常接近.Matsui 等[11]采用Zienkiewicz 等[10]的方法分析多個邊坡的穩(wěn)定性,并把該方法正式命名為強(qiáng)度折減法.Griffiths 等[12]詳細(xì)論述了如何將強(qiáng)度折減法與理想彈塑性有限元法結(jié)合分析邊坡的穩(wěn)定性,并通過大量算例分析及與極限平衡法結(jié)果比較,說明有限元強(qiáng)度折減法的有效性.有限元強(qiáng)度折減法中安全系數(shù)大小與采用的巖土材料屈服準(zhǔn)則密切相關(guān)[13].傳統(tǒng)的極限平衡法采用摩爾?庫侖準(zhǔn)則,但是由于摩爾?庫侖準(zhǔn)則的屈服面是截面為不規(guī)則的六角形的角錐體表面,存在尖頂和菱角,給數(shù)值計算帶來困難.一些商業(yè)軟件采用廣義米賽斯準(zhǔn)則,或摩爾?庫倫的外接圓屈服面.這些與傳統(tǒng)的摩爾?庫倫準(zhǔn)則的計算結(jié)果存在較大差異[13].徐干成和鄭穎人[14]提出的等面積圓屈服面用于強(qiáng)度折減法計算邊坡的穩(wěn)定安全系數(shù),與Spencer 法計算的誤差在5%左右[13],表明采用合適的本構(gòu)模型,強(qiáng)度折減法不但可以找到最危險滑動面位置,還能定量地計算穩(wěn)定安全系數(shù).強(qiáng)度折減法中一般采用計算不收斂作為邊坡臨界的判據(jù)[15],由于有限元彈塑性迭代計算的復(fù)雜性,計算結(jié)果與軟件本身的彈塑性迭代算法關(guān)系較大.筆者發(fā)現(xiàn)利用ABAQUS(版本6.14-1)軟件采用摩爾?庫倫模型對復(fù)雜的邊坡進(jìn)行計算時,甚至得不到收斂的安全系數(shù)結(jié)果; 對于簡單邊坡淺表層滑動機(jī)制,也得不到合理的安全系數(shù),如找不到坡角為30?單一砂土(無凝聚力)邊坡正確的塑性極限狀態(tài).雖然強(qiáng)度折減法是邊坡穩(wěn)定分析的強(qiáng)有力的工具,在國內(nèi)土質(zhì)邊坡和巖質(zhì)邊坡的穩(wěn)定性中也有較多的應(yīng)用,但強(qiáng)度折減法所得的安全系數(shù)至今仍尚未被國內(nèi)相關(guān)規(guī)范所采用.另外,該方法對最危險滑動面以外的其他潛在滑動面的安全系數(shù)也難以分析.

極限分析用于邊坡穩(wěn)定分析,是 Drucker 和Prager (1952) 首先完成的.他們用滿足平衡條件且在任意一點都不違背屈服準(zhǔn)則的假想應(yīng)力場,確定了破壞載荷的下限.上限是通過一個與流動法則相容的速度場,利用虛功原理,即外力的虛功率等于或大于內(nèi)部的能量耗損率[16].陳惠發(fā)介紹了利用上限法求解直線、圓弧和對數(shù)螺線滑動面的邊坡穩(wěn)定極限分析方法[16].極限分析法與條分法(包括垂直條分法和斜條分法) 相結(jié)合,通過求解虛功方程,可獲得給定組合剛體滑動機(jī)制的邊坡穩(wěn)定安全系數(shù).并通過優(yōu)化方法,求得臨界的剛體組合滑動機(jī)制和最小上限安全系數(shù)[17].條分法由于其固有的弱點,難以適應(yīng)復(fù)雜的地質(zhì)情況和復(fù)雜的滲流與變形耦合作用問題.極限分析與剛體有限元法結(jié)合的方法,也得到了很大發(fā)展.系統(tǒng)的塑性變形內(nèi)能耗散僅發(fā)生在單元間的界面上,故此類方法的計算網(wǎng)格,需要適應(yīng)滑動面的形狀,在優(yōu)化過程中需要進(jìn)行自適應(yīng)網(wǎng)格重構(gòu)[18],大大增加了程序算法的難度.極限分析法還存在難以考慮復(fù)雜的滲流應(yīng)力耦合相互作用的局限性.

有限元方法求解邊坡的應(yīng)力場,可以考慮任意復(fù)雜的地質(zhì)組成和滲流變形耦合情況.將有限元應(yīng)力結(jié)果與極限平衡法分析的概念結(jié)合起來的方法,稱為改進(jìn)的極限平衡法[19].Kulhawy[20]1969 年提出的改進(jìn)極限平衡法用一個獨立的程序求出應(yīng)力場后,通過插值的方法獲得滑動面上點的應(yīng)力狀態(tài),再用滑面上抗剪強(qiáng)度的代數(shù)和與剪切力的代數(shù)和之商計算穩(wěn)定安全系數(shù),能求解復(fù)雜地質(zhì)和復(fù)雜滲流變形耦合情況下邊坡的穩(wěn)定性.這一方法對于圓弧滑動,由于滑弧上每一點切向力對于圓心的力臂均相同,實際上這抗滑力與滑動力在滑面上的代數(shù)和之比等于抗滑力矩與滑動力矩之比,因而是合理的.而對于非圓弧和非直線滑動面,力的代數(shù)和概念不清,存在問題.對于任意滑動面,劉艷章等[21]將滑面上法向力和抗剪強(qiáng)度的合力定義為抗滑力矢,并將整個滑面上抗滑力矢的和方向定義為主滑動方向,提出采用給定滑動面上總抗滑力矢在主滑動方向上投影的代數(shù)和與總下滑力矢在此方向上投影代數(shù)和的比值定義安全系數(shù),進(jìn)行了若干工程應(yīng)用研究.然而,該方法對于圓弧滑動機(jī)制,其安全系數(shù)與常規(guī)極限平衡方法顯然不一致.對于實際發(fā)生的圓弧與折線組合式滑動或多段折線組合式滑動,其主滑面方向的物理含義也不明確,因而,該方法所求出的安全系數(shù)的物理意義不明確.

本文假定邊坡的滑動機(jī)制為組合剛體滑動,基于有限元應(yīng)力場與運(yùn)動許可的組合剛體速度場,依據(jù)虛功原理,提出了邊坡穩(wěn)定安全系數(shù)計算的新方法,探討了該方法所基于的有效應(yīng)力場是否需要滿足靜力許可的問題,并在典型算例中與其他極限平衡方法分析結(jié)果進(jìn)行了對比.

1 邊坡穩(wěn)定分析的虛功率法理論

極限分析法通過計算靜力許可的應(yīng)力場或機(jī)動容許的速度場,根據(jù)虛功原理來計算極限載荷.上限定理可表述為:與機(jī)動容許速度場對應(yīng)的外載荷不小于真實的極限載荷.物體所能承受的最大外載荷,與極限狀態(tài)下虛速度作用下物體內(nèi)部的虛功率成正比.虛應(yīng)變能可以作為表征物體承載力的單一指標(biāo).下面根據(jù)邊坡在機(jī)動容許的虛速度場下的虛功率來定義和計算安全系數(shù).

1.1 安全系數(shù)的定義與計算

非單一圓弧或非單一直線滑動面,包括折線、非圓弧曲線、線段與圓弧的組合等,統(tǒng)稱為組合式滑動面.如果將滑體作為剛體,則不存在沿著組合式滑動面的機(jī)動容許的速度場的單一剛體,因此,組合式滑動面上的邊坡體是剛體的組合.圖1 是與斜條分極限分析法[17]中相同的滑動機(jī)制,邊坡被幾個速度間斷面分割成若干剛體.速度間斷面可以為圓弧、直線或?qū)?shù)螺線.圖中ABCDEF折線稱為底滑面或主滑面,A,B,C,D,E,F為主滑點,2 個端點A,F稱為坡面主滑點,其他主滑點稱為內(nèi)部主滑點,BB1和CC1等稱為內(nèi)部錯動面或次滑面,B1和C1等稱為次滑點.下文所稱的滑動面,包括主滑面和次滑面.

圖1 Fig.1

邊坡滑動時,滑面上存在相對速度.滑動體總的功率等于滑動面上的功率和,可以表示為

其中,P為功率,τ 和σn分別為滑動面上的剪應(yīng)力和正應(yīng)力;v為滑動面上的相對速度; ? 為運(yùn)動速度與間斷面之間的夾角,因剪脹而引起,又稱之為剪脹角,不同滑動面上的剪脹角可以不相同;Γ 為所有速度間斷面.

邊坡的滑動臨界狀態(tài)通過抗剪強(qiáng)度的折減達(dá)到.剛體極限平衡法和有限元強(qiáng)度折減法都用臨界狀態(tài)時的強(qiáng)度折減系數(shù)的倒數(shù)作為安全系數(shù).假定邊坡速度間斷面上的正應(yīng)力不受土體抗剪強(qiáng)度折減的影響,臨界狀態(tài)時滑動面上的虛功率與強(qiáng)度未折減時依據(jù)有限元計算所得應(yīng)力場的間斷面上的總虛功率相等,設(shè)安全系數(shù)為F,則

其中τf為滑面上的抗剪強(qiáng)度,采用摩爾庫倫強(qiáng)度準(zhǔn)則

其中c為凝聚力,φ 為內(nèi)摩擦角.

由式(2) 可以得到給定組合剛體邊坡的滑動安全系數(shù)F的顯式表達(dá)為

上式即安全系數(shù)等于滑動面上總的抗剪強(qiáng)度功率與總的剪應(yīng)力功率之比.分子和分母均是標(biāo)量求和,物理含義明確.

虛功率法基于機(jī)動容許的速度場和有限元應(yīng)力場,所得的安全系數(shù)是給定滑動機(jī)制的上限解.邊坡穩(wěn)定分析中邊坡的應(yīng)力計算是一個獨立的步驟,因而可以不受限制地采用獨立的有限元分析法獲得任意復(fù)雜地質(zhì)情況與耦合情況下的有效應(yīng)力應(yīng)力場.

速度間斷面的位置可以通過分步優(yōu)化的方法來確定,以獲得最危險滑動機(jī)制.求邊坡給定底滑面的穩(wěn)定安全系數(shù),是內(nèi)部錯動面的位置經(jīng)過優(yōu)化后所得的最小安全系數(shù).而整個邊坡的安全系數(shù),則是所有可能的主滑面中,安全系數(shù)最小的那個.

1.2 滑動面上剪切力與抗剪強(qiáng)度的計算

二維有限元網(wǎng)格將滑動面分割為若干線段(或弧段),可以首先計算滑動面上與有限元網(wǎng)格的交點,找出這些線段及其所在的單元.依據(jù)有限元計算所得的每個單元的高斯點應(yīng)力,通過等參逆變換插值法求出每個單元的節(jié)點應(yīng)力.滑動面上任意一點的應(yīng)力可以通過單元的節(jié)點應(yīng)力插值計算,即

式中,σi j為所求點處的應(yīng)力張量,下標(biāo)i,j表示笛卡爾坐標(biāo)系空間坐標(biāo)軸;NJ為計算點關(guān)于節(jié)點J的形函數(shù);為單元內(nèi)節(jié)點編號為J處的應(yīng)力張量.

滑面積分點上的應(yīng)力向量Ti可由下式計算

式中,nj表示滑面積分點處法線方向在整體坐標(biāo)中的方向向量.法向應(yīng)力和剪應(yīng)力分別為

其中l(wèi)i為切向方向向量.滑動面高斯積分點的抗剪強(qiáng)度采用摩爾庫倫公式(3)計算.

1.3 邊坡組合滑動機(jī)制與速度場

在有限元計算獲得邊坡的應(yīng)力場以后,通過構(gòu)造組合剛體機(jī)動容許的速度場,就可以計算得到給定滑動機(jī)制下的安全系數(shù).如圖1 所示的組合剛體滑動機(jī)制,主滑面是A,B,C,D,E,F點構(gòu)成的5 段折線.主滑面也可以由折線、圓弧(剪脹角為0)和對數(shù)螺線(剪脹角不為0)構(gòu)成.圖中AB,BC,CD,DE,EF,BB1,CC1,DD1,EE1為速度間斷面.速度間斷面上的運(yùn)動方向與間斷面成? 角(剪脹角),不考慮剪脹時?=0.設(shè)第一個滑入段的速度為單位速度1,根據(jù)速度三角形可以求出各間斷面上的速度.如圖1(b)所示,速度向量構(gòu)成的速度三角形ABC中,滑面上的相對速度的方向與x軸正方向的夾角等于滑動面的方向角加上(或減去)一個剪脹角.可以計算出速度三角形的3 個內(nèi)角,∠A,∠B和∠C.根據(jù)三角形邊角關(guān)系定理,可得

1.4 安全系數(shù)的計算與優(yōu)化方法簡介

根據(jù)滑面上的力和速度,求得每一條速度間斷面上的滑動力功率和抗滑力功率后,依據(jù)式(4)即可計算得到給定滑動機(jī)制的邊坡穩(wěn)定安全系數(shù).

尋找臨界滑動面的優(yōu)化方法是多變量優(yōu)化問題.虛功率法中的優(yōu)化變量包括所有主滑點和次滑點的位置.多變量優(yōu)化中變量越多,優(yōu)化難度越大.采取分步優(yōu)化的方法可降低難度.分成3 層嵌套:最內(nèi)層是主滑點不動次滑點優(yōu)化; 次內(nèi)層是坡面主滑點不動內(nèi)部主滑點優(yōu)化.

最內(nèi)層次滑面位置優(yōu)化也是一個多變量優(yōu)化問題.對次滑點位置可以采用多輪單變量順序優(yōu)化,將多變量優(yōu)化問題,化簡為多次單變量優(yōu)化問題.次內(nèi)層內(nèi)部主滑點位置的優(yōu)化也是多變量優(yōu)化問題.可以按照文獻(xiàn)[19]提出的下降差分方法和收斂準(zhǔn)則,采用梯度法準(zhǔn)確找到極小值位置.對于位于軟弱結(jié)構(gòu)面上的內(nèi)部點,可以規(guī)定只沿著結(jié)構(gòu)面移動,對于其他情況,可以沿著次滑面方向移動.最外層坡面主滑點的位置沿著坡面線移動,其優(yōu)化是2 變量優(yōu)化問題,與內(nèi)部主滑點優(yōu)化方法相同.

2 邊坡穩(wěn)定安全系數(shù)計算檢驗

作者開發(fā)的LinkFEA-Slope 程序?qū)⒁陨侠碚摲椒{入其中,程序具有基于有限元有效應(yīng)力計算給定圓弧或折線組合式滑動面安全系數(shù)的功能,并采用梯度法優(yōu)化尋找最危險滑動面.

1987 年Donald 教授和Giam 博士主持了澳大利亞計算機(jī)協(xié)會(ACADS)對澳大利亞所使用的邊坡穩(wěn)定分析程序進(jìn)行了調(diào)查研究.共設(shè)計了5 個考題總計10 個小題,向120 個單位發(fā)出測試邀請,28 個單位發(fā)回了計算答案.Donald 教授也使用自編的GWEDWEM 和EMU 程序分析了考題,并請以色列Baker 教授用SSA 程序、中國的陳祖煜教授用STAB 程序和加拿大的Fredlund 教授提供裁判答案.最終綜合各單位提供的答案、裁判答案和Donald 教授的計算結(jié)果給出了推薦答案.其中的第3 題(EX3)的第2 小題和第4 題(EX4)是測試非圓弧滑動的,用來檢驗本文理論方法的合理性.

2.1 算例介紹

EX3 考題的邊坡幾何剖面和有限元網(wǎng)格見圖2.材料參數(shù)見表1.地下水位位于軟弱夾層(2#土) 底部,地下水位以上孔壓為零.第1 問為圓弧或折線滑動面的最小安全系數(shù)(EX3-1),第二問為計算指定的折線滑動面ABCD(控制點坐標(biāo)如表2) 的安全系數(shù)(EX3-2).

圖2 EX3 邊坡幾何剖面及指定的折線滑動面與有限元網(wǎng)格Fig.2 Slope geometry,the specified polygonal sliding surface and the finite element mesh in EX3

表1 EX3 材料參數(shù)Table 1 Material parameters in EX3

表2 EX3 指定折線滑動面控制點Table 2 Control points of the specified polyline sliding surface

考題EX4 的邊坡如圖3 所示,坡面上作用有豎向載荷、坡內(nèi)有軟弱夾層和地下水.邊坡材料參數(shù)見表3,浸潤線描述見表4.要求確定臨界滑動面和計算相應(yīng)最小安全系數(shù).

圖3 EX4 邊坡剖面及載荷情況及有限元網(wǎng)格Fig.3 Slope profile,load condition and the finite element mesh in EX4

表3 EX4 材料參數(shù)Table 3 Material parameters in EX4

表4 EX4 浸潤線坐標(biāo)Table 4 Phreatic line position in EX4

2.2 應(yīng)力場計算條件

應(yīng)力場采用ABAQUS(版本6.14)計算.采用四邊形孔壓?應(yīng)力耦合單元(CPE4P).首先約束所有節(jié)點的位移,求出穩(wěn)定滲流場,然后將負(fù)孔壓節(jié)點的孔壓修改為0 后,節(jié)點孔壓作為已知條件,采用線彈性本構(gòu)模型求解有效應(yīng)力場.

EX3 的滲流場比較簡單,是浸潤線位于弱夾層底部,該部位以下孔隙水壓力為靜水壓力,其上為0孔隙水壓力.EX4 原題只提供了浸潤線,而有限元應(yīng)力計算必須要有節(jié)點孔隙水壓力.有限元滲流計算中夾層的滲透系數(shù)取1#土滲透系數(shù)的50 倍,右邊界取水面高程38.4 m 的等水頭邊界條件,左邊界及臺地地表取27.75 m 水頭邊界,斜坡取透水邊界條件,所得邊坡滲流場浸潤線與等孔壓線如圖4 所示,圖中可見浸潤線位置與原題所給有一定差異.

圖4 EX4 邊坡浸潤線與等孔壓線Fig.4 The phreatic line and pressure contour in EX4

2.3 安全系數(shù)結(jié)果與文獻(xiàn)結(jié)果的比較

應(yīng)力場基于線彈性參數(shù)計算結(jié)果的EX3-2 與EX4 的安全系數(shù)結(jié)果與文獻(xiàn)中給出的裁判答案列于表5.

表5 算例安全系數(shù)結(jié)果比較Table 5 Comparison of safety coefficient results of examples

文獻(xiàn)中29 個提交的答案(同一研究者依據(jù)不同極限平衡方法提交多個答案)EX3-2 安全系數(shù)平均值為1.292,標(biāo)準(zhǔn)差(均方差)為0.064,推薦答案為1.340.本文EX3-2 中給定折線滑動面的安全系數(shù)1.360,比推薦答案大1.5%,介于裁判答案之間.EX4 文獻(xiàn)中19個提交答案的平均安全系數(shù)為0.817,標(biāo)準(zhǔn)差0.223,推薦的裁判答案為0.78[23].本文EX3 給定折線主滑動面、優(yōu)化次滑面位置的安全系數(shù)如圖5 所示;EX4的初始和優(yōu)化后的折線主滑面折點坐標(biāo)如表6.EX4折線滑動面如圖6,B點次滑面方位角127.0?,C點次滑面方位角75.4?.安全系數(shù)0.703,介于裁判答案之間,比推薦答案小9.9%.2 個含有軟弱夾層的算例計算結(jié)果都介于裁判答案之間,與推薦答案接近,表明虛功率法對于折線組合式滑動方案的研究結(jié)果合理.

圖5 EX3 滑動面的位置Fig.5 Position of the sliding surface in EX3

表6 EX4 滑動面初始點與最優(yōu)點(括號中)坐標(biāo)Table 6 Coordinates of initial and optimum(in parentheses)points on the sliding surface

圖6 EX4 優(yōu)化后的滑動面位置Fig.6 The optical sliding surface position in EX4

采用ABAQUS(版本6.14)軟件中的強(qiáng)度折減法,本構(gòu)模型采用M-C 模型,單元類型采用2 次三角形單元,取不同的網(wǎng)格密度分別計算EX3 與EX4 的穩(wěn)定安全系數(shù),發(fā)現(xiàn)塑性區(qū)圖所顯示的滑動機(jī)制基本一致,而網(wǎng)格越密,安全系數(shù)越小.0.1 m 邊長網(wǎng)格模型計算所得的塑性區(qū)與穩(wěn)定安全系數(shù)如圖7 所示,安全系數(shù)分別為1.261 和0.785(網(wǎng)格邊長1.0 m 的模型安全系數(shù)分別為1.285 和0.810).虛功率法應(yīng)力計算的網(wǎng)格尺度大于1.0 m.對比來看,虛功率法的速度間斷面位置與有限元強(qiáng)度折減法所得的結(jié)果基本一致.而安全系數(shù)結(jié)果算例EX3 基本一致,EX4 小約10%.

圖7 強(qiáng)度折減法計算塑性區(qū)與安全系數(shù)Fig.7 The plastic zone and safety factor calculated by strength reduction method

2.4 穩(wěn)定安全計算方法討論

依據(jù)極限分析的下限定理,有限元強(qiáng)度折減法的應(yīng)力場是靜力許可的,極限狀態(tài)所對應(yīng)的外載荷小于真實的極限承載力,其安全系數(shù)結(jié)果理論上是一個下限值.由于非線性迭代的收斂性和網(wǎng)格密度的影響,計算所得的極限狀態(tài)時的應(yīng)力場并不能處處滿足靜力許可.且由于塑性應(yīng)變是局部化的,因而一般需要采用細(xì)密網(wǎng)格和二次單元才能得到比較準(zhǔn)確的結(jié)果(一次單元對網(wǎng)格密度要求更高).網(wǎng)格越細(xì)密,程序的非線性功能越好,計算結(jié)果越接近這個下限值.其缺點之一是要求程序具有強(qiáng)大的非線性計算能力,這種能力往往在地質(zhì)和滲流變形耦合情況下難以達(dá)到;其二是需要細(xì)密的網(wǎng)格,對于復(fù)雜的工程問題計算規(guī)模往往過大而成本過高,甚至計算任務(wù)難以完成.筆者還未見采用商業(yè)軟件用強(qiáng)度折減法求滲流與變形雙向耦合情況下的穩(wěn)定安全系數(shù)的工程案例介紹.

極限平衡法對于同一土條內(nèi)包含多種土層的情況誤差較大,且難以考慮復(fù)雜的滲流性狀,尤其是難以考慮滲流應(yīng)力強(qiáng)耦合相互作用的情況.而基于有限元應(yīng)力計算的改進(jìn)極限平衡法對于復(fù)雜滑動面安全系數(shù)計算又存在明顯的理論缺陷.

虛功率法應(yīng)力場計算獨立于后續(xù)的穩(wěn)定性分析,采用有限元法完成,復(fù)雜的地質(zhì)和滲流應(yīng)力耦合作用由有限元完成,因而不受限制.應(yīng)力計算一次完成,因而速度間斷面位置優(yōu)化的效率也遠(yuǎn)遠(yuǎn)高于極限平衡法.復(fù)雜滑動機(jī)制又繼承了極限分析上限法的優(yōu)點,因而對于解決地質(zhì)條件和滲流應(yīng)力耦合作用問題有顯著的優(yōu)點.與有限元強(qiáng)度折減法相比,一方面應(yīng)力計算時對于材料非線性迭代計算的性能要求較低,另一方面對于網(wǎng)格密度要求較低(不需要依靠計算形成塑性應(yīng)變連通區(qū)得到穩(wěn)定安全系數(shù)).計算的效率較高.

虛功率法安全系數(shù)計算公式來源于極限分析的上限法,其結(jié)果是一個上限解.結(jié)果與實際的接近程度,取決于假定的滑動機(jī)制的合理性及滑動機(jī)制的優(yōu)化效果.當(dāng)然,對于邊坡穩(wěn)定安全系數(shù)接近或小于1.0的情況,線彈性模型計算所得的應(yīng)力場僅滿足靜力平衡,違反靜力許可的區(qū)域較大,安全系數(shù)計算的結(jié)果與實際值差異相對較大.對于實際安全系數(shù)小于1.0的情況,是不存在靜力許可的應(yīng)力場的.而對于安全裕度不大的邊坡,通過采用彈塑性模型,或?qū)`反強(qiáng)度準(zhǔn)則的單元,進(jìn)行應(yīng)力遷移等非線性迭代計算,減小違反靜力許可的區(qū)域、降低單元區(qū)域應(yīng)力超過破壞面的程度,則安全系數(shù)結(jié)果會更接近于實際.

2.5 應(yīng)力場結(jié)果差異對安全系數(shù)的影響

邊坡的有限元應(yīng)力場計算結(jié)果,受土的本構(gòu)模型及其參數(shù)影響.基于以下條件的計算獲得EX3 算例4 個應(yīng)力場結(jié)果:(1)線彈性模型和如表1 中的參數(shù)(A); (2)線彈性模型且2#土也取表1 中1#土的參數(shù)(B);(3)摩爾庫倫彈塑性模型和表1 中的參數(shù)(C);(4)摩爾庫倫模型且強(qiáng)度除以1.25 折減(D).

基于不同的應(yīng)力場,計算了EX3 給定滑動面位置的安全系數(shù)、也計算了最危險圓弧滑動面的安全系數(shù),滑動面位置如圖5 所示.對于給定的主滑面,次滑面位置優(yōu)化后方位角如表7.BB1方位角最大差異1.5?,CC1方位角最大差異4.2?,優(yōu)化后不同應(yīng)力場的次滑面位置存在一定差異但差異不太大.

表7 EX3-2 應(yīng)力條件與最優(yōu)次滑面方位角Table 7 Stress condition and azimuth angle of optimal subslip surfaces

給定折線主滑面的安全系數(shù)和圓弧滑動面最小安全系數(shù)結(jié)果比較分別如表8 和表9.折線滑動面安全系數(shù)介于1.333～1.367 之間,最大差異0.034,相對于A 條件的幅度為?2.0%.本例對于圓弧滑動,A 應(yīng)力條件下,最危險圓弧位置如圖6,滑入、滑出和圓心點的坐標(biāo)分別為(72.45,40)、(42.54,27.75)、(50.35,51.32).不同應(yīng)力條件下圓弧滑動安全系數(shù)最大差異幅度為0.9%.可見本例中本構(gòu)模型和參數(shù)造成的應(yīng)力場差異對安全系數(shù)的影響很小.

表8 應(yīng)力條件與折線滑動面安全系數(shù)Table 8 Stress condition and safety factor of polyline sliding surface

表9 應(yīng)力條件與臨界圓弧滑動面安全系數(shù)Table 9 Stress condition and safety factor of polyline sliding surface

彈性模型計算得到的應(yīng)力場均滿足靜力平衡,但不一定滿足靜力許可; 彈塑性模型計算所得的應(yīng)力場,既滿足靜力平衡,也基本滿足靜力許可.應(yīng)力條件不同,虛功率法計算得到的安全系數(shù)結(jié)果存在差異.根據(jù)最小值原理,一個機(jī)構(gòu)的極限承載力,不低于任意一個同時滿足靜力平衡和靜力許可的應(yīng)力極限狀態(tài)所對應(yīng)的外載荷.因此對于應(yīng)力場既滿足靜力平衡,又滿足靜力許可的邊坡,無疑是穩(wěn)定的,其安全系數(shù)不低于1.即使是穩(wěn)定的邊坡,當(dāng)安全裕度不高時,無論是采用什么樣的本構(gòu)模型,有限元應(yīng)力計算時常常存在高斯點應(yīng)力超過其抗拉和抗剪強(qiáng)度的情況,而要通過彈塑性迭代或應(yīng)力遷移迭代使應(yīng)力結(jié)果滿足強(qiáng)度準(zhǔn)則,則往往是十分困難的.邊坡穩(wěn)定安全系數(shù)小于1.0 的情況下,根本就不存在靜力許可的應(yīng)力場.要進(jìn)行強(qiáng)度折減(提高)后才能獲取靜力許可的應(yīng)力場.因此,對于安全系數(shù)計算所依據(jù)的應(yīng)力場,僅要求靜力平衡而不要求靜力許可一般還是合理的.采用虛功率法,依據(jù)靜力平衡的應(yīng)力場計算所得的穩(wěn)定安全系數(shù)是邊坡穩(wěn)定性較好的度量.

3 結(jié)論和展望

本文提出了邊坡穩(wěn)定安全性計算的虛功率法理論和實施方法,在2 個典型折線滑動邊坡算例中進(jìn)行了檢驗.并探討了應(yīng)力場靜力平衡與靜力許可條件對安全系數(shù)結(jié)果的影響.得出以下結(jié)論:

(1)虛功率法用邊坡組合剛體滑動機(jī)動許可的速度間斷面上抗滑力功率與滑動力功率之比計算安全系數(shù),解決了改進(jìn)的極限平衡法中對于復(fù)雜滑動面安全系數(shù)計算公式物理含義不清晰的問題;

(2) 采用分步優(yōu)化方法對組合滑動機(jī)制的速度間斷面的位置優(yōu)化,能較容易準(zhǔn)確確定最危險滑動機(jī)制;

(3)虛功率法的應(yīng)力計算獨立于安全系數(shù)的計算與優(yōu)化過程,因而能考慮復(fù)雜的地質(zhì)和滲流與應(yīng)力耦合作用情況;

(4)采用靜力許可的應(yīng)力場,安全系數(shù)結(jié)果更接近實際,靜力平衡的應(yīng)力場,其穩(wěn)定安全性計算結(jié)果,也是邊坡穩(wěn)定性的一個不錯的度量.

虛功率法是為適應(yīng)水電站庫岸堆積體高邊坡和復(fù)雜軟弱水工構(gòu)筑物地基的穩(wěn)定性分析提出的新的穩(wěn)定分析方法.包含該方法的分析軟件LinkFEASlope 軟件已經(jīng)用于瀾滄江如美和金沙江拉哇2 個重大水電站工程的庫岸堆積體邊坡、碎裂卸荷巖體邊坡或深厚軟弱覆蓋層高圍堰的穩(wěn)定性分析,正在用于長江沿岸某物流碼頭軟土地基上超高堆場基礎(chǔ)穩(wěn)定性分析.采用虛功率法,基于有限元應(yīng)力場結(jié)果分析圓弧與非圓弧滑動機(jī)制的邊坡、地基穩(wěn)定安全系數(shù),得到了工程設(shè)計部門和業(yè)主的歡迎與鼓勵.目前作者研究團(tuán)隊一方面加緊改進(jìn)復(fù)雜的滲流與變形耦合相互作用的計算程序LinkFEA-Stress,意圖應(yīng)力計算不但包含滲流與應(yīng)力的強(qiáng)耦合作用,也包含拉應(yīng)力遷移和塑性滑移在內(nèi)的材料非線性模擬,力圖使應(yīng)力場更接近靜力許可; 另一方面加緊開發(fā)包括圓弧與折線和對數(shù)螺線管等復(fù)雜滑弧組合式底滑面滑動機(jī)制的安全系數(shù)計算與優(yōu)化,將在工程應(yīng)用中不斷改進(jìn)與發(fā)展.