非均質邊坡穩定性有限元強度折減影響因素分析

馬 偉

(山東農業大學勘察設計研究院，山東 泰安 271000)

1 引 言

土石壩、堤防、水庫庫岸等水利工程中存在著大量的邊坡問題。邊坡一旦失穩破壞，將造成重大的人員傷亡和財產損失[1]。因此，邊坡穩定性分析一直是土石構筑物工程中的重要研究領域。近年來，將強度折減法和有限元結合應用于邊坡穩定性分析已受到了越來越多的關注。強度折減有限元法克服了傳統極限平衡法的諸多缺點，能考慮土體的非線性本構關系，滿足力的平衡條件，模擬復雜的邊界條件和多種材料；不僅能直觀地得出邊坡的應力、應變、特征點位移、塑性區分布等參數值及相應的云圖分布，而且能直接得到一個物理意義明確的定量評價邊坡穩定性的安全系數。

但就目前而言，在強度折減有限元計算中，單元類型、網格大小的選擇，初始地應力考慮與否，邊坡失穩評價標準的選取等都會對最后的計算結果造成影響；巖土體的黏聚力、內摩擦角、重度、彈性模量、泊松比以及剪脹角等材料參數同樣會對安全系數的計算產生影響。

針對均質邊坡的研究較多。萬少石等[2]討論了單元類型、網格大小對安全系數計算精度的影響，各種失穩判據的一致性、土體材料參數和坡角對邊坡破壞模式的影響，得出了一些有益的結論。王棟等[3]通過比較計算認為單元類型應選四邊形二階單元。徐文杰等[4]討論了初始地應力對有限元強度折減計算結果的影響。薛雷等[5]對非均質邊坡穩定性中強度折減的范圍進行了研究，認為應對整個區域進行折減。張魯渝等[6]討論了坡高、坡角、黏聚力、摩擦角對強度折減有限元法計算精度的影響，并給出了提高計算精度的具體措施。邊坡的失穩判據主要有3類：?以數值計算迭代不收斂為失穩判據；?以特征部位的位移突變為失穩標志；?以是否形成連續的塑性貫通區為失穩判據。對于上述3類失穩判據，不同學者分別從各種失穩判據的差異性和內在統一性進行了論述。

目前，將有限元強度折減應用于均質邊坡的例子較多，但很少涉及非均質邊坡，本文擬研究單元類型和網格大小、初始地應力、失穩判據等因素對非均質邊坡強度折減的影響。

2 強度折減法基本原理

20世紀70年代，Zienkiewics首先提出了強度折減的概念，后被許多學者廣泛采用，并成功地應用于實際工程中。強度折減法就是在理想彈塑性有限元計算中，通過不斷折減土體的抗剪強度參數c、φ，直至邊坡達到臨界破壞狀態，此時的強度折減系數就是邊坡的穩定安全系數Fs，定義為土體實際的抗剪強度與折減后極限狀況下抗剪強度的比值，通過式(1)、式(2)進行計算：

c′=c/F

(1)

φ′=arctan(tanφ/F)

(2)

式中c——土體實際的黏聚力，kPa；

φ——土體實際的摩擦角，(°)；

c'——土體折減后的黏聚力，kPa；

φ'——土體折減后的的摩擦角，(°)；

F——強度折減系數。

本文采用有限元分析程序進行有限元強度折減，在非均質邊坡中對整個區域進行折減[7]，本構模型采用Mohr-Coulomb屈服準則和非關聯流動法則。

3 強度折減有限元影響因素分析

3.1 單元類型和網格大小的影響

采用強度折減有限元法分析非均質邊坡穩定性時，合理地選取網格的單元類型和大小是保證強度折減計算精度的重要條件之一。本文以非均質邊坡[8]為例進行強度折減有限元分析，并與簡化Bishop法所得結果比較，探討網格單元類型和大小對非均質邊坡強度折減有限元計算精度的影響。

非均質邊坡的剖面幾何尺寸及材料分區見圖1，材料特性見表1。

圖1 非均質邊坡的剖面幾何尺寸及材料分區 (單位：m)

表1 非均質邊坡的材料特性

為研究不同單元類型和網格大小對計算精度的影響，分別選取了四邊形和三角形的一階、二階單元，網格大小為0.5～8m，每隔0.5m變化，進行強度折減有限元分析，計算出了在不同單元類型和網格大小下的安全系數，網格劃分采用自由劃分技術。為方便討論，以有限元數值計算不收斂為失穩判據。強度折減所得安全系數和簡化Bishop法的計算結果見圖2。該算例的簡化Bishop法所得結果為0.43[8]。

圖2 單元類型和網格大小對計算精度的影響

由圖2可知，不管采用何種單元類型，計算結果均隨網格的稀疏而增大。與簡化Bishop法所得結果相比，四邊形一階單元網格大小為0.5～3m時，誤差在5%以內，當網格大小為1.5～2.5m時，誤差僅在2%以內；四邊形二階單元網格大小為0.5～5.5m時，誤差在5%以內，當網格大小為1.5～4.5m時，誤差僅在2%以內；采用三角形一階單元只有網格大小在0.5m時，誤差在5%以內，網格變大，則得不到較好的計算精度；三角形二階單元網格大小為1～9m時，誤差在5%以內，當網格大小為2～5.5m時，誤差僅在2%以內。由分析可知，采用三角形一階單元時，計算精度得不到保證，一般不宜采用；采用四邊形一階、二階單元時，只要網格大小合理均能保證計算結果的精度；三角形二階單元對網格大小的要求不高，即使網格較大，仍能獲得較理想的結果。由此可知，采用有限元強度折減對非均質邊坡進行穩定性分析是可行的，對于四邊形單元，當相對網格大小(網格大小與坡高的比值)為0.075～0.125；三角形單元，相對網格大小為0.1～0.275時，所得安全系數與Bishop法所得結果較為接近。

在非均質邊坡強度折減有限元中，并不是網格劃分得越小，所得結果越準確。二階單元時，網格劃分過小，不但耗費的機時很長，而且結果有時還不可用。比較一階單元和二階單元，二階單元所得安全系數隨網格大小變化較小，結果的穩定性好，建議采用二階單元分析非均質邊坡的穩定性。

3.2 初始地應力的影響

為了研究初始地應力對非均質邊坡穩定性的影響，分別采用四邊形一階、二階單元和三角形二階單元，在考慮初始地應力情況下進行強度折減。根據計算結果，以數值計算不收斂為判據，不考慮初始地應力和考慮初始地應力的計算結果見表2。

由表2可知，當采用四邊形單元時，考慮初始地應力時所得結果偏小，采用三角形二階單元時，是否考慮初始地應力對計算結果沒有影響，對于非均質邊坡的強度折減有限元分析，初始地應力對所得結果有一定的影響，條件允許時要考慮初始地應力的影響。

表2 安全系數計算結果

3.3 邊坡失穩判據分析

如前所述，在邊坡穩定性分析中應用強度折減有限元主要有3種判據，失穩判據的選取一直以來都有較大爭議。本例采用三角形二階單元4m網格進行邊坡穩定分析，將坡頂點水平位移與強度折減系數的關系曲線繪于圖3中，圖4給出了不同強度折減系數下的塑性區分布。

圖3 坡頂點水平位移與強度折減系數的變化曲線

圖4 塑性區隨強度折減系數的變化

由圖3可知，若以坡頂點水平位移突變為失穩判據，可得安全系數為Fs=0.423，對比圖4中塑性區分布，當強度折減系數F=0.423時，邊坡塑性區開始出現貫通現象，表征邊坡正好處于極限狀態。由此可得，采用特征點位移突變和塑性區是否貫通所得安全系數是一致的。若以數值計算不收斂為失穩判據，安全系數為Fs=0.426，所得結果與前兩種判據的計算結果相比，數值偏大。分析邊坡失穩破壞的過程，隨強度折減的進行，邊坡內土體開始出現塑性區并不斷發展，擴大直至貫通，伴隨著滑動體位移發生突變，而后進入無限制的塑性流動，繼而無法承受荷載。強度折減有限元數值計算不收斂與所選取的屈服準則、所采用的計算模型、單元類型和網格大小以及設定的迭代步數等因素有關，數值計算不收斂并不一定意味著邊坡達到極限破壞狀態，或者可能邊坡早已失穩。因此，雖然在本例中以坡頂點位移突變和以數值計算不收斂所得到的計算結果十分接近，但這與選取的單元類型和網格大小有關，并不具有普遍性。采用特征點位移突變結合塑性區貫通作為失穩判據具有明確的物理意義，更能真實反映邊坡的穩定性。

4 算例驗證

選取澳大利亞計算機應用協會(ACADS)對邊坡穩定分析程序調查時使用的一個考核題[9]。邊坡的剖面幾何尺寸及材料分區見圖5，材料特性見表3。

圖5 算例邊坡的幾何尺寸及材料分區

表3 算例邊坡的材料特性

有限元強度折減時，采用三角形二階單元，2.5m網格，自由劃分技術。將考慮初始地應力下所得安全系數與ACADS計算結果見表4。

表4 計算結果匯總表

由表4可知，Bishop法所得結果較大，Janbu法所得結果較小，都偏離裁判答案。通過有限元強度折減以坡頂點位移突變為失穩判據所得結果1.394與裁判答案1.39幾乎沒有誤差，可見本文所推薦的單元類型和網格大小以及失穩判據具有一定的合理性，對于簡單的非均質邊坡是適用的。但對于更為復雜的非均質邊坡，結論的普遍性有待進一步研究。

5 結 語

本文將強度折減法與有限元結合應用于非均質邊坡的穩定性分析。討論了單元類型和網格大小、初始地應力、失穩判據對非均質邊坡有限元強度折減計算結果的影響，得到了以下結論:

a.在非均質邊坡強度折減有限元中，網格采用三角形一階單元時計算精度差，四邊形一階、二階單元在合理網格大小范圍內滿足精度要求，三角形二階單元對網格大小要求不高，精度高。建議采用三角形二階單元分析非均質邊坡，參考的合理相對網格大小為0.1～0.275。

b.初始地應力對計算結果有一定的影響，考慮初始地應力所得結果較不考慮時偏小。若條件允許，建議在非均質邊坡穩定性分析中考慮初始地應力的影響。

c.分析表明，采用特征點位移突變結合塑性區貫通為失穩判據具有明確的物理意義，能真實反映邊坡的穩定性。