王少華 許小明 蔡 瑾
(湖北三峽職業技術學院,湖北宜昌 443000)
邊坡穩定性分析理論經過多年的研究發展已經逐步成熟,極限平衡法主要表現在其力學模型簡單,無法反應模型內部應力應變情況以及破壞機理,且在安全系數計算方面操作繁瑣;有限元法雖然克服了極限平衡法以上各方面的不足,也能計算各種復雜邊界的邊坡,但不能和安全系數建立直接聯系而給出一個量化的評價標準,在一定程度上限制了其應用范圍。本文主要以重力加載法確定的潛在軟弱面為指定滑弧,計算土坡極限平衡狀態下的安全系數。
ANSYS程序是一個功能強大、靈活的設計分析優化分析軟件包,可以浮動運行于從PC,NT工作站直至巨型機的各類計算機及操作系統中,與眾多先進CAD軟件共享數據接口,數據文件在其所有的產品系列和工作平臺上均兼容。
Slide是一款評價巖質或土質邊坡安全系數或者失效概率的二維極限平衡程序,滑面兼具圓弧和非圓弧形式,可指定已知滑面或者驅動程序使之自動搜索滑面。
本文在研究過程中均采用ANSYS,Slide進行分析。
本文假定被研究模型為理想彈塑性材料,其遵循關聯流動法則,更符合巖土體的實際情況。
大變形有限元重度增加法與極限平衡綜合的基本思路分為三步[1-3]:第一步,通過大變形有限元法重度增加法確定邊坡的最危險滑動面及坡體內應力應變情況;第二步,按極限平衡條分法原理將滑坡體分條,分析條塊的受力情況,并推導安全系數公式以及編制計算程序;第三步,根據最危險滑動面位置和形狀擬合成滑動面方程,將滑動面方程以及巖土體材料輸入計算程序,精確求解安全系數。研究方法用結構示意圖可以簡單的表示,如圖1所示。
邊坡滑動面以及潛在滑動面[6]作為邊坡穩定分析的一個重要環節,如從傳統極限平衡法的角度來看,確定邊坡臨界滑動面時,搜索范圍的指定具有任意性,并且是基于一定假設的基礎之上,不一定能夠搜索到真正的危險滑移面。
根據本文思想,如通過數值模擬來確定滑動面位置所在,則失穩準則的選擇在整個分析過程中顯得尤為重要。判定邊坡失穩的理論方法目前概括起來只有三種:1)有限元數值迭代計算不收斂;2)特征部位特征點位移突變;3)等效塑性應變從坡腳至坡頂的貫通。由于塑性區貫通是破壞的必要非充分條件,那么,危險滑移面必定發生在塑性區貫通面上,本文嘗試通過有限元重力加載法找到該塑性區貫通面。

圖1 研究方法結構示意圖
具體問題具體分析,以下用模型實例來探討在重力增加的過程中,塑性區的擴展(漸變)過程。算例基本資料:淹鍋沙壩某均質土坡,邊坡模型以及物理力學參數如表1所示。

表1 物理力學參數表
通過大變形有限元重度增加法在有限元軟件ANSYS中模擬邊坡破壞過程:隨著重力加速度的增大,坡內塑性區不斷重分布并不斷增大,最大塑性應變逐漸形成塑性應變等值脊線密集區,并不斷向上發展直至貫通,形成明顯的滑動帶,如圖2中MX標志處。在重度增加到1.82g時,塑性區已基本貫通,但未形成滑移面,在重度增加至1.84g時,塑性區域內形成了貫通的滑移面,但不能判斷是否已經發生了大變形,而在重度增加到1.86g時,滑動面上部開始閉合,說明大變形已經發生,故在1.82~1.84之間的任意一個數值都可能作為超載儲備安全系數,如果人為觀察確定塑性區剛剛貫通時的數值這個操作是難以實現的,所以單以塑性區貫通為標志來確定安全系數值是不合理的,從而也間接證明了重度增加法的塑性區貫通判據只能作為邊坡失穩的必要非充分條件,與前人研究結果一致。形成貫通的滑移面即為潛在滑動面,可以作為極限平衡計算中的指定滑帶。
對于圓弧滑動面,簡化Bishop法是具有足夠精度的工程實用計算方法;對于任意形狀滑動面,只有嚴格按條分法才能得到合理的安全系數,而且條間力函數的變化對安全系數值影響也很小。鑒于此,本文參考極限平衡法中簡化Bishop法來推導安全系數計算公式不失一般性。

圖2 不同重度幅值下邊坡內塑性區發展及分布圖
下面針對邊坡來進行計算分析,驗證本文方法的可靠性。
工程背景:淹鍋沙壩某土坡,坡高35 m,坡長90 m,坡角45°,假三維有限元模型共有6 813個單元,9 026個節點,模型邊界采用固定約束,坡體材料選用M-C模型,材料參數見表1。
由于土體所發生的位移遠遠小于邊坡幾何尺寸,這是經典有限元常用的小變形假設[9]。在有限元軟件計算中打開大變形分析,選擇大變形關聯流動法則來模擬分析,加載子步100,最大循環50次,收斂容差0.05。通過編程計算得出各自結果見圖3~圖6。

圖3 強度折減FS=1.243塑性區貫通圖

圖4 極限平衡FS=1.116失穩滑動面圖
圖7中1號滑弧為有限元重度增加法確定的滑面,2號滑弧為有限元強度折減法確定的滑面,3號滑弧為極限平衡法確定的滑面。由于邊坡的最危險滑動面應穿過塑性區上最大塑性應變的峰值點,由此借助AutoCAD尋找滑弧圓心(-4.2,63.85)半徑63.987,以及滑弧剪出口坐標(0,0),(81.341,50)擬合最危險滑動面方程,輸入已編制的程序,計算安全系數。從以上三條滑動面位置對比發現最危險滑面位置基本一致,只是由于重度增加法與強度折減法得到的塑性區大小不同,主要是施加荷載的方式不同引起的坡體內部網格單元達到屈服狀態的先后順序不一而導致的。

圖5 原重度增加FS=1.432塑性區貫通圖

圖6 原重度增加法塑性應變等值線圖

圖7 CAD描圖三條滑動面位置對比

表2 不同方法計算安全系數的對比
從表2安全系數上對比分析,綜合法與傳統極限平衡簡化Bishop法的結果最為接近,與強度折減法的結果相比次之,而與原重度增加法的結果相差最大。宏觀評價以上各種方法對安全系數精度的影響,極限平衡法是建立在對最危險滑面形狀的一系列假定,并視巖土體為剛體的前提之下,通過復雜的內置計算程序搜索最小安全系數對應的滑面,其本身與實際情況就有一定的出入。而有限元法自身的精度影響和收斂容差的設置,會出現安全系數計算誤差,而綜合法兼顧兩者之長,取長避短。表3驗證了運用綜合法分析的結果較為理想,并得出原重度增加法的分析結果誤差較大。

表3 最危險滑面確定的難易度比較
1)本文方法不但使滑動面的確定簡單明了,且避免了以往傳統極限平衡法中對滑動面位置的假設、搜索和安全系數試算等一系列復雜的過程,同時避免了有限元軟件中命令流、計算程序等復雜的調試過程以及人為設置帶來的偏差。2)有限元重度增加法和極限平衡法的綜合法在分析求解安全系數方面思路簡潔清晰,且能夠建立安全系數與穩定分析的直接聯系,具有一定的參考價值。3)有限元重度增加法在計算安全系數方面存在較大的誤差,本文分析重度增加法安全系數誤差主要來源于重度在滑面方向的投影分力,并因此對重度增加法安全系數計算公式加以改進,減小了誤差。4)依然從重度增加法安全系數計算式建立抗滑力增量與下滑力之間的不等式,從而確定該方法的適用范圍。
[1]謝永利.大變形固結理論及其有限元法[M].北京:人民交通出版社,1998.
[2]任青文.非線性有限單元法[M].南京:河海大學出版社,2003.
[3]唐春安,唐烈先.巖土破裂過程分析RFPA離心加載法[J].巖土工程學報,2007,29(1):31-32.
[4]鄧建輝.CAD技術在邊坡穩定性分析中的應用[D].武漢:中國科學院武漢巖土力學研究所,1992.
[5]林 鵬,卓家壽.邊坡穩定性廣角度分析的等值面法[J].水利學報,2000(8):75-79.
[6]王成華,夏緒勇,李廣信.基于應力場的土坡臨界滑動面的螞蟻算法搜索技術[J].巖石力學與工程學報,2003,22(5):813-819.
[7]傅薈璇.MATLAB應用設計[M].北京:機械工業出版社,2010.
[8]Hibbitt H D,Marcal P V,Rice J R.A finite element formulation for problems of large strain and large displacement[J].International Journal of Solids and Structure,1970(6):1069-1086.
[9]孫廣忠.工程地質理論與實踐[M].北京:科學出版社,1997.