5 剛體單元上限法

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

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

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

5.1 非均質土坡旋轉破壞機構

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

圖5-1中三段對數螺旋線滑裂面具有相同的轉動中心O，并且相鄰兩層土交界處，不同對數螺旋線的半徑保持不變。每一層土對應的對數螺旋線滑弧，其函數關系可以用內摩擦角表示為

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

當坡高H和坡度β確定后，通過坡趾的滑裂面Α0Α1Α2C，完全可以由具有共同旋轉中心O的分段對數螺旋線確定；而旋轉中心O可以由x和y的距離確定，或由θ0和θ3單獨確定。

圖5-1 非均質土坡旋轉破壞機制及速度矢量

5.2 層狀土質邊坡外力功率計算

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

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

圖5-2 滑體三角形單元網格化簡圖

OOei與x軸正方向夾角為

設旋轉破壞機構角速度為ω，則編號為i的三角形單元應變速度為

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

同理，可以計算得到所有單元在應變速度場v所做的外功率為

圖5-3 剛體單元計算外力功率

5.3 對數螺旋線組合滑面內能耗散計算

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

同理，可以計算得到第二層土對數螺旋線Α1Α2段和第三層土對數螺旋線Α2C段的內能耗散分別為：

根據數學歸納法，可以得到第i層土的內能耗散計算公式

邊坡滑體包括多層土時，則組合對數螺旋線總內能耗散為：

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

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

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

式中：

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

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

以某一兩層土組成的非均質土坡為例，坡高8m，坡比1∶1.3，從坡頂到坡腳土體分兩層，分別記為第①層和第②層，厚度各4m，土層參數見表5-1。通過采用剛體極限平衡法（對應圓弧滑面）和剛體單元上限法（對應分段對數螺旋線滑面）兩種不同方法搜索最危險滑移面，得到圖5-4中顯示的結果。

表5-1 雙層土質邊坡土層參數

圖5-4 雙層土質邊坡剖面

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

表5-2 不同方法的穩(wěn)定系數計算結果

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

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

表5-3 地震敏感性分析經驗函數

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