吳夢(mèng)喜 楊家修 湛正剛

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

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

??(中國(guó)電建集團(tuán)貴陽(yáng)水利水電勘測(cè)設(shè)計(jì)研究院,貴陽(yáng) 550081)

引言

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

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

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

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

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

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

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

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

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

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

圖1 Fig.1

邊坡滑動(dòng)時(shí),滑面上存在相對(duì)速度.滑動(dòng)體總的功率等于滑動(dòng)面上的功率和,可以表示為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 邊坡組合滑動(dòng)機(jī)制與速度場(chǎng)

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

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

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

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

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

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

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

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

2.1 算例介紹

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

圖2 EX3 邊坡幾何剖面及指定的折線(xiàn)滑動(dòng)面與有限元網(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 指定折線(xiàn)滑動(dòng)面控制點(diǎn)Table 2 Control points of the specified polyline sliding surface

考題EX4 的邊坡如圖3 所示,坡面上作用有豎向載荷、坡內(nèi)有軟弱夾層和地下水.邊坡材料參數(shù)見(jiàn)表3,浸潤(rùn)線(xiàn)描述見(jiàn)表4.要求確定臨界滑動(dòng)面和計(jì)算相應(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 浸潤(rùn)線(xiàn)坐標(biāo)Table 4 Phreatic line position in EX4

2.2 應(yīng)力場(chǎng)計(jì)算條件

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

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

圖4 EX4 邊坡浸潤(rùn)線(xiàn)與等孔壓線(xiàn)Fig.4 The phreatic line and pressure contour in EX4

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

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

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

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

圖5 EX3 滑動(dòng)面的位置Fig.5 Position of the sliding surface in EX3

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論和展望

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

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

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

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

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

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