抗滑樁設計推力計算的有限元方法

雷 達，蔣關魯，林展展，張樹明，王智猛

(1.西南交通大學土木工程學院，四川成都 610031;2.西南交通大學 道路與鐵道工程鐵道部開放實驗室，四川成都 610031;3.中鐵二院工程集團有限責任公司，四川 成都 610031)

抗滑樁加固滑坡已經成為滑坡災害治理的主流工程措施，它充分利用抗滑樁和樁周巖土體的交互作用，平衡滑坡推力，同時分擔部分下滑推力到基巖錨固段。目前，利用傳遞系數法和相應修正公式進行抗滑樁設計的方法，因其簡單實用的力學原理和成熟的工程應用而成為主流方法。潘家錚［1］提出，先計算達到設計安全系數時滑坡斷面的下滑推力，再計算自然條件下的下滑推力，差值為抗滑樁設計推力。由于忽略巖土體自身穩定能力，計算結果偏于保守。賀建清等［2］認為抗滑樁在受荷的同時對樁前土體或圍巖產生反力，假定樁前抗力傳遞系數為a，設計推力必須減去樁前抗力，但取值條件有待商榷。和海芳等［3］認為抗滑樁前后滑體擁有相等的安全系數，提出以滑坡剪出口附近滑塊推力是否收斂于0為判別條件，確定樁前剩余下滑力，以抗滑樁前后推力差值作為設計標準。王培勇等［4］分別從上往下計算滑坡推力，從下往上計算抗滑力，以兩者差值作為抗滑樁設計推力。

以上計算遵守剛體極限平衡假定，滑面滿足靜力平衡，但不滿足力矩平衡。只能進行近似分析，不能真實表現土的彈塑性應力—應變關系，抗滑樁設計推力計算存在一定誤差。程建軍等［5］研究表明有限元法的滑坡推力計算結果相對合理，應該推廣應用。本文結合成蘭鐵路右所屯滑坡，利用有限元軟件FLAC3D模擬滑坡地應力場、位移場、邊界條件等，進行非線性分析，并與傳遞系數法進行比較。

1 有限元法計算原理

主要用到兩個本構計算模型:各向同性彈性模型、摩爾—庫倫模型［6-7］。低應力水平時使用各向同性彈性模型，認為土體近似產生各向同性彈性變形，根據胡克定律有

式中:K，G分別為體積模量、剪切模量，δij為克羅內克δ函數，Δσij為應力張量增量，Δεij，Δεkk為應變張量增量，應力計算沒有殘余變形分量。使用摩爾—庫倫模型模擬滑坡塑性流動破壞時，彈性變形產生的廣義應力增量仍然用胡克定律進行計算。

滑坡破壞是由降雨、地震等工況引起的滑面強度參數的降低，彈性模型不涉及。利用摩爾—庫倫模型得到的剩余下滑力應扣除彈性解，即:F剩余下滑力=

2 有限元法算例

2.1 工點概況

右所屯工點位于成蘭鐵路 D3K247+788—D3K247+846段，滑體軸向長100 m，寬40～60 m，滑體厚3～15 m，為一中型堆積層滑坡，線路從滑坡前緣通過。滑體為粉質黏土，滑帶為細角礫土，基巖為以弱風化為主的炭質板巖夾砂巖。根據前期地質調查結果，滑體和滑帶之間存在潛在滑動面，在降雨、地震等惡劣工況下，將沿滑面發生整體滑動破壞。為簡化計算分析，滑體粉質黏土換為細角礫土，滑坡土層物理參數如表1所示。

表1 滑坡土層物理參數

2.2 建立數值計算模型

監測斷面劃分如圖1所示，共8個，監測斷面受力為網格節點合力，滑體和滑帶分界面附近存在10 cm厚細角礫土潛在滑動面。根據右所屯工點滑坡主軸斷面圖，建立單位寬度的有限元計算模型，共有單元體2 542個、節點5 284個。在滑坡模型的X方向兩側邊界和Y方向，用fix設置固定節點速度和位移的邊界條件，底部邊界同理為X，Y，Z三向的速度位移固定約束。

圖1 監測截面示意(單位:m)

2.3 剩余下滑力分步計算

計算前，移除監測截面沿滑坡剪出口方向的巖土體。以x=64 m監測截面為例，移除監測截面至左側邊界之間的巖土體，并對x=64 m巖土體邊界施加固定節點位移和速度的約束，按照以下三步計算:①利用彈塑性求解法生成初始地應力，鑒于黏聚力為0，對滑面摩擦角進行折減，在摩爾—庫倫本構模型下數值計算得到受力穩定后的截面合力F1，運行時步5萬步。②利用彈性本構模型生成初始地應力，得到計算截面的彈性合力F3。③由于滑坡失穩的根本原因是摩爾—庫倫本構模型中，潛在滑動面的摩擦角變小產生下滑力［8］，而彈性本構模型只涉及體積模量、剪切模量、密度，故在推力中扣除彈性解部分，F1和F3的差值F4為滑面摩擦角折減時監測斷面的剩余下滑力。

傳遞系數法計算公式見式(9)和(10)，為便于和傳遞系數法的計算結果進行比較，分為4種工況分別計算剩余下滑力，即潛在滑動面內摩擦角分別折減為30°，25°，20°，10°。

式中:Ei為第i個條塊末端的滑坡推力，kN/m;ks為安全系數(視工程重要性、整治難易程度);Wi為第i個條塊滑體重力，kN/m;αi－1為第 i－1個條塊所在滑動面的傾角，(°);αi，φi，ci，ui，li分別為第 i個條塊滑動面的傾角、內摩擦角、單位黏聚力、單位孔隙水壓力和長度。

2.4 結果分析

圖2為有限元法和傳遞系數法剩余下滑力對比，為便于分析，剩余下滑力均取X方向(水平方向)分力。傳遞系數法下滑推力變化的線性特征明顯，有限元法曲線變化緩和平穩，能較好地模擬推力的非線性變化。隨著滑動面土層摩擦角減小，滑動面摩阻力變小，滑體內部彈塑性作用效果減弱，推力曲線線性特征逐漸明顯。當滑動面摩擦角為10°時，兩種算法推力曲線的線性特征非常顯著。傳遞系數法在滑面摩擦角大幅折減時，抗滑段下滑推力沒有收斂減小趨勢，一直傳遞到最后一個土塊;有限元法則能清晰分辨下滑段和抗滑段。當滑動面摩擦角為30°時，傳遞系數法推力曲線在距離坡腳26 m時已經收斂為0，這與實際工況不符。有限元法則能顯示滑坡推力先增加后減小的規律，沒有截面推力突變為0的情況。由于兩種算法的理論基礎不同，傳遞系數法的推力峰值大于有限元法，曲線變化規律一致，下滑力峰值都出現在X方向88 m處。

圖2 剩余下滑力對比

圖3 剩余下滑力增量對比

圖3為滑動面摩擦角變化時的剩余下滑力增量對比。當滑動面摩擦角變化區間為30°～25°，25°～20°時，兩種方法的剩余下滑力增量變化趨勢一致，且傳遞系數法增量均大于有限元法，說明傳遞系數法對剩余下滑力增量有放大效應，這種方法對工程設計偏于安全。隨著滑動面摩擦角折減加劇，理論上25°～20°區間的剩余下滑力增量應大于30°～25°時的，但傳遞系數法下滑段的5個監測截面不滿足此規律。當滑動面摩擦角變化區間為20°～10°時，抗滑力急劇減少導致滑體產生劇烈位移，兩種方法的增量變化在下滑段相近。但在抗滑段表現出不同的趨勢，傳遞系數法的下滑力增量不斷增加，沒有收斂趨勢。

2.5 設計推力計算

自然狀態時儲備的抗滑力使滑體處于穩定狀態。隨著滑面摩擦角不斷減小，當滑坡處于極限破壞平衡狀態時，臨界摩擦角出現，抗滑力已發揮全部作用，監測截面的剩余下滑力和抗滑力平衡。有限元法利用FISH程序語言編寫LOOP循環語句［9］，對滑面摩擦角進行折減試算。以程序內置力不平衡比率＞1×10－5為終止標準，得到滑面臨界摩擦角為27.63°，安全系數為1.338，計算顯示滑坡上部表層發生溜塌。傳遞系數法按照式(9)和式(10)試算，直到滑坡剪出口附近滑塊推力近似為0，得到滑面臨界摩擦角為24.85°，安全系數為1.512。

以滑面摩擦角20°時為基準，比較抗滑樁設計推力。有限元法按照剩余下滑力的計算步驟，得到臨界極限狀態和基準條件時的剩余下滑力，兩者差值即為抗滑樁設計推力。

2.6 結果分析

圖4顯示兩種算法的設計推力曲線有一致的變化規律，但傳遞系數法各點設計推力均大于有限元法，這種差異在抗滑段達到最大。

圖4 設計推力曲線

3 結論

1)和傳遞系數法相比，利用有限元法模擬塑性流動和破壞，曲線變化相對平緩，能縮小剩余下滑力的計算誤差，不會出現抗滑段剩余下滑力計算不收斂或提前收斂于0的情況。

2)有限元法的設計推力變化趨勢與傳遞系數法一致，且小于傳遞系數法的計算結果。

［1］潘家錚.建筑物的抗滑穩定和滑坡分析［M］.北京:水利出版社，1980:155-184.

［2］賀建清，張家生，梅松華.彈性抗滑樁設計中幾個問題的探討［J］.巖石力學與工程學報，1999，18(5):600-602.

［3］和海芳，祁生文，伍法權，等.抗滑樁設計推力計算方法研究［J］.工程地質學報，2008，16(5):694-698.

［4］王培勇，劉元雪，冉仕平，等.基于傳遞系數法確定支護樁設計推力的新方法［J］.四川大學學報(工程科學版)，2010，42(6):73-78.

［5］程建軍，廖小平，王浩，等.滑坡推力計算方法的對比研究與應用［J］.水文地質工程地質，2008(1):44-47.

［6］陳育民，徐鼎平.FLAC/FLAC3D基礎與工程實例［M］.北京:中國水利水電出版社，2009.

［7］Itasca Consulting Group，Inc..Fast Language Analysis of Continua in 3 Dimensions，Version 3.0，User's Manual［M］.Itasca Consulting Group，Inc.，2005.

［8］DAWSON E M，ROTH W H，DRCSCHER A.Slope Stability Analysis by Strength Reduction［J］.Geotechnique，1999，49(6):835-840.

［9］彭文斌.FLAC3D實用教程［M］.北京:機械工業出版社，2008.