5 剛體單元上限法

土坡穩(wěn)定性分析是巖土工程經(jīng)典問(wèn)題之一，在鐵路、公路、水利水電、建筑基坑及港口碼頭工程中應(yīng)用非常廣泛。土坡穩(wěn)定性計(jì)算常用的方法是剛體極限平衡法，如Fellenius（1927）提出的圓弧滑動(dòng)分析方法（瑞典圓弧法），以及改進(jìn)的畢肖普法（1955）、斯賓塞法（1967）、簡(jiǎn)布法（1973）、摩根斯坦-普賴斯法（1965）以及薩爾瑪法（1973）等。盡管剛體極限平衡法以其簡(jiǎn)單實(shí)用而得到了廣泛的應(yīng)用，但因此方法本身存在一定的局限性，直接導(dǎo)致所取得的解并非總是正確，故其不能解決土體破壞的運(yùn)動(dòng)學(xué)問(wèn)題。相反，基于土塑性力學(xué)的極限分析上限法，卻可以選擇機(jī)動(dòng)許可的速度場(chǎng)來(lái)尋求土坡的破壞機(jī)制（陳祖煜，2002），并且這一方法已被廣泛認(rèn)可和采用。

通過(guò)采用極限分析上限法，陳惠發(fā)（1975）最早提出了解決均質(zhì)土坡臨界高度的方法，并采用類似土壓力系數(shù)的穩(wěn)定數(shù)Ns來(lái)進(jìn)行土坡穩(wěn)定性分析。但這一方法沒(méi)能推廣到非均質(zhì)土坡問(wèn)題，也沒(méi)能給出通用的穩(wěn)定系數(shù)計(jì)算公式。Karal（1977）后來(lái)對(duì)陳惠發(fā)的方法進(jìn)行了改進(jìn)，得到了旋轉(zhuǎn)破壞機(jī)制直接分析評(píng)價(jià)均質(zhì)土坡的穩(wěn)定系數(shù)計(jì)算公式。Michalowski（1995）則根據(jù)極限分析上限法提出了類似條分法的平動(dòng)破壞機(jī)制，并得到了均質(zhì)土坡計(jì)算穩(wěn)定數(shù)Ns的公式。Kumar 和 Samui（2006）還根據(jù)非相關(guān)聯(lián)流動(dòng)法則，采用虛功率方程得出了土坡穩(wěn)定系數(shù)計(jì)算公式。

在土質(zhì)邊坡極限分析上限法領(lǐng)域，存在兩個(gè)明顯的問(wèn)題：其一是不能考慮成層土的非均質(zhì)特性，其二是沒(méi)有通用的穩(wěn)定系數(shù)計(jì)算公式。正是基于目前存在的問(wèn)題，筆者提出了將極限分析上限法和有限單元法前處理功能相融合的思想，并通過(guò)強(qiáng)度折減法和虛功率方程，實(shí)現(xiàn)了求解非均質(zhì)土坡穩(wěn)定系數(shù)的數(shù)值計(jì)算方法。

5.1 非均質(zhì)土坡旋轉(zhuǎn)破壞機(jī)構(gòu)

均質(zhì)土坡破壞時(shí)的機(jī)動(dòng)許可速度場(chǎng)（破壞機(jī)構(gòu)），根據(jù)極限分析法可知其可能產(chǎn)生的滑面為一對(duì)數(shù)螺旋線。對(duì)于二維三層土組成的非均質(zhì)土坡，每一層土可以通過(guò)層厚hi、黏聚力ci、內(nèi)摩擦角φi和容重γi四個(gè)參數(shù)進(jìn)行定義（i=1，2，3），其可能滑裂面為通過(guò)坡趾的不同對(duì)數(shù)螺旋線的組合（見圖5-1）。

圖5-1中三段對(duì)數(shù)螺旋線滑裂面具有相同的轉(zhuǎn)動(dòng)中心O，并且相鄰兩層土交界處，不同對(duì)數(shù)螺旋線的半徑保持不變。每一層土對(duì)應(yīng)的對(duì)數(shù)螺旋線滑弧，其函數(shù)關(guān)系可以用內(nèi)摩擦角表示為

式中：r0，r1，r2，r3和θ0，θ1，θ2，θ3等符號(hào)意義參見圖5-1，且有

當(dāng)坡高H和坡度β確定后，通過(guò)坡趾的滑裂面Α0Α1Α2C，完全可以由具有共同旋轉(zhuǎn)中心O的分段對(duì)數(shù)螺旋線確定；而旋轉(zhuǎn)中心O可以由x和y的距離確定，或由θ0和θ3單獨(dú)確定。

圖5-1 非均質(zhì)土坡旋轉(zhuǎn)破壞機(jī)制及速度矢量

5.2 層狀土質(zhì)邊坡外力功率計(jì)算

陳惠發(fā)（1975）在其著作《極限分析與土體塑性》中，推導(dǎo)了均質(zhì)土坡外力功率計(jì)算公式。但對(duì)于非均質(zhì)土坡，常規(guī)的解析方法計(jì)算外力功率難以實(shí)現(xiàn)。筆者根據(jù)有限元軟件ΑNSYS的前處理功能——網(wǎng)格化，通過(guò)三節(jié)點(diǎn)三角形單元網(wǎng)格化非均質(zhì)土體，可以分別計(jì)算每一個(gè)單元的外力功率。

圖5-2所示為滑體經(jīng)由三角形單元網(wǎng)格化后的效果。從網(wǎng)格化后的滑體中取出編號(hào)為i的三角形單元，如圖5-3（a）所示，該單元節(jié)點(diǎn)編號(hào)按逆時(shí)針依次為j，k，l，節(jié)點(diǎn)坐標(biāo)分別為（xi，j，yi，j），（xi，k，yi，k），（xi，l，yi，l）。三角形單元（element-i，記為ei）面積計(jì)算公式為

圖5-2 滑體三角形單元網(wǎng)格化簡(jiǎn)圖

OOei與x軸正方向夾角為

設(shè)旋轉(zhuǎn)破壞機(jī)構(gòu)角速度為ω，則編號(hào)為i的三角形單元應(yīng)變速度為

編號(hào)為i的三角形單元自重為Wei=γei△ei，所受到的水平地震作用力為khWei，如圖5-3（b）所示。外力作用在應(yīng)變速度vei所做的外功率為

同理，可以計(jì)算得到所有單元在應(yīng)變速度場(chǎng)v所做的外功率為

圖5-3 剛體單元計(jì)算外力功率

5.3 對(duì)數(shù)螺旋線組合滑面內(nèi)能耗散計(jì)算

根據(jù)極限分析上限法對(duì)土質(zhì)邊坡穩(wěn)定分析的假定，滑體將被視為剛體，而滑裂面被視為速度間斷面。因此在滑體內(nèi)將不存在內(nèi)能的耗散，換言之，旋轉(zhuǎn)機(jī)構(gòu)的內(nèi)能耗散集中于組合對(duì)數(shù)螺旋線上。以圖5-1所示的三層土質(zhì)邊坡為例，首先計(jì)算第一層土對(duì)數(shù)螺旋線Α0Α1段的內(nèi)能耗散，取微分長(zhǎng)度rdθ/cosφ1，此長(zhǎng)度內(nèi)的黏聚力為c1rdθ/cosφ1，則黏聚力沿滑裂面的總內(nèi)能耗散積分表達(dá)式為

同理，可以計(jì)算得到第二層土對(duì)數(shù)螺旋線Α1Α2段和第三層土對(duì)數(shù)螺旋線Α2C段的內(nèi)能耗散分別為：

根據(jù)數(shù)學(xué)歸納法，可以得到第i層土的內(nèi)能耗散計(jì)算公式

邊坡滑體包括多層土?xí)r，則組合對(duì)數(shù)螺旋線總內(nèi)能耗散為：

5.4 剛體單元法穩(wěn)定系數(shù)上限解

在極限狀態(tài)下，多層土組成的非均質(zhì)邊坡虛功率方程為

將式（5-11）、式（5-16）和式（5-18）代入上式，可得

式中：

式（5-21）即為求解非均質(zhì)土坡穩(wěn)定性穩(wěn)定系數(shù)的剛體單元上限法計(jì)算公式。由于穩(wěn)定系數(shù)K是以隱式出現(xiàn)的，因此需要通過(guò)迭代的算法求解。在計(jì)算時(shí)，先假設(shè)一個(gè)初始的穩(wěn)定系數(shù)K，應(yīng)用迭代計(jì)算直至前后兩次穩(wěn)定系數(shù)K值滿足所要求的計(jì)算精度為止。一般來(lái)說(shuō)對(duì)于滿足條件的情況，迭代都是收斂的。

采用剛體單元上限法計(jì)算非均質(zhì)土坡穩(wěn)定系數(shù)的步驟為：① 通過(guò)ΑNSYS有限元軟件對(duì)非均質(zhì)邊坡建模，并進(jìn)行三角形單元剖分；② 輸出節(jié)點(diǎn)文件（包括節(jié)點(diǎn)號(hào)和節(jié)點(diǎn)坐標(biāo)）和單元文件（包括單元編號(hào)和單元對(duì)應(yīng)節(jié)點(diǎn)號(hào)）；③ 編程計(jì)算所有單元在重力及地震作用力下產(chǎn)生的外力功率；④ 采用強(qiáng)度折減法，計(jì)算每層土對(duì)應(yīng)的對(duì)數(shù)螺旋線內(nèi)能耗散；⑤ 通過(guò)式（5-21）多次迭代獲得剛體單元上限法穩(wěn)定系數(shù)最終解。

以某一兩層土組成的非均質(zhì)土坡為例，坡高8m，坡比1∶1.3，從坡頂?shù)狡履_土體分兩層，分別記為第①層和第②層，厚度各4m，土層參數(shù)見表5-1。通過(guò)采用剛體極限平衡法（對(duì)應(yīng)圓弧滑面）和剛體單元上限法（對(duì)應(yīng)分段對(duì)數(shù)螺旋線滑面）兩種不同方法搜索最危險(xiǎn)滑移面，得到圖5-4中顯示的結(jié)果。

表5-1 雙層土質(zhì)邊坡土層參數(shù)

圖5-4 雙層土質(zhì)邊坡剖面

剛體極限平衡法分別選用瑞典條分法、畢肖普法和簡(jiǎn)布法，在自重作用下，穩(wěn)定系數(shù)分別為0.96，1.00和1.00（表5-2）。在相同條件下，剛體單元上限法得到的穩(wěn)定系數(shù)為1.01。從計(jì)算結(jié)果來(lái)分析，兩者均可以得到邊坡處于臨界破壞狀態(tài)的結(jié)論，但文中方法略比剛體極限平衡法結(jié)果要大。同樣，考慮Ⅶ度地震烈度、Ⅷ度地震烈度和Ⅸ度地震烈度，計(jì)算結(jié)果均顯示邊坡將產(chǎn)生破壞，但穩(wěn)定系數(shù)也呈現(xiàn)比剛體極限平衡法略大。

表5-2 不同方法的穩(wěn)定系數(shù)計(jì)算結(jié)果

通過(guò)四種方法得到的穩(wěn)定系數(shù)與地震烈度的敏感性分析來(lái)看（見圖5-5），文中方法曲線斜率要小于剛體極限平衡法。采用最小二乘法對(duì)以上數(shù)據(jù)進(jìn)行函數(shù)擬合，得到表5-3的穩(wěn)定系數(shù)——地震烈度的經(jīng)驗(yàn)函數(shù)關(guān)系。綜合分析，可以得到三點(diǎn)結(jié)論：①剛體單元上限法穩(wěn)定系數(shù)對(duì)地震烈度的敏感性要小于其他三種方法；②隨著地震烈度的增加，四種方法得到的穩(wěn)定系數(shù)對(duì)地震烈度的敏感性均有所提高；③畢肖普法地震烈度敏感性最大。分析計(jì)算結(jié)果顯示，文中方法和畢肖普法可以代表穩(wěn)定系數(shù)的上限和下限。因此，通過(guò)多種方法計(jì)算穩(wěn)定系數(shù)時(shí)，在界定了上限和下限后，合理選用穩(wěn)定系數(shù)對(duì)邊坡工程非常重要。

圖5-5 地震烈度敏感性分析

表5-3 地震敏感性分析經(jīng)驗(yàn)函數(shù)

剛體單元上限法的主要局限性在于它從上限逼近真實(shí)解，因此，會(huì)過(guò)高地估計(jì)穩(wěn)定系數(shù)，這在工程上使用是偏于不安全的。畢肖普法得到的穩(wěn)定系數(shù)相對(duì)要小，這在工程上使用是偏于安全的。綜合這幾種方法的計(jì)算結(jié)果，穩(wěn)定系數(shù)的真實(shí)解應(yīng)該介于剛體單元上限法和畢肖普法兩者之間。