考慮土體各向異性的抗滑樁邊坡穩定性分析

寧 帥 莊 妍 談云志 王東星

（1.東南大學 蘇州聯合研究生院, 南京 211189；2.東南大學 土木工程學院, 南京 211189；3.三峽大學三峽庫區地質災害教育部重點實驗室, 湖北 宜昌 443002；4.武漢大學 土木建筑工程學院, 武漢 430072）

關于土拱效應的研究一直受到許多學者的關注,1943年,Terzaghi(太沙基)[1]率先驗證了土拱效應的真實存在,ITO等[2]通過塑性流動理論和塑性變形理論推導出土拱效應下樁的側向力公式,Richard等[3]針對小主應力的作用軌跡,提出了土拱形狀是一條懸鏈線的觀點；Kellogg等[4]將土拱按大致的形狀分為半球形、三角形、圓頂形、拋物線形等.由于土拱形狀的不確定性和土體參數的復雜性,導致關于土拱效應的理論計算方法十分有限,并且大多數的計算方法都是以一定的假設為前提條件.通過數值模擬的方法,學者們在土拱效應的研究上取得了許多進展.Ladanyi等[5-10]利用理論與數值模擬相結合的方法,對土拱效應展開了更為詳細的研究；董捷等[11]運用數值模擬對三維土拱效應進行了研究并將土拱分為水平、豎向和臨空3類.目前比較常見的研究是通過數值模擬得到的數據與理論計算得出的結果進行對比驗證,對現有的計算方法進行更準確的修正.同時,部分學者將對土拱的研究結果應用于室內試驗與實際工程中,為工程的設計和安全提供了更加充分地保障[12-15].在抗滑樁加固邊坡的穩定性研究中,已有許多學者對抗滑樁后土體的土拱效應進行了詳細的研究[16-21],通過這些研究可以發現,土體幾乎都被視為各向同性,并且簡化了數值模擬試驗的模型和加載過程,因此無法充分體現出真實情況下抗滑樁加固邊坡的復雜性.在抗滑樁邊坡的數值模擬試驗中,極少有學者能夠在有限元軟件中考慮到各向異性來模擬實際工況.然而,就筆者調查研究所知,如果在試驗過程中將土體視為各向同性,計算結果所得的安全系數會明顯高于土體為各向異性時的安全系數,并且差距會隨著土體的參數和模型的尺寸而改變.因此,在計算抗滑樁邊坡穩定性時,需要將土體的各向異性條件考慮在內.王棟等[22]利用ABAQUS和強度折減法實現了各向異性條件下邊坡安全系數的數值分析,但是其模擬的情況為未加固的普通邊坡,極少有學者在有限元軟件中充分考慮各向異性土體和水平土拱效應來研究抗滑樁加固邊坡的穩定性.本文針對以上問題建立了抗滑樁加固邊坡的力學模型,通過對ABAQUS子程序的編寫及調用,根據迭代計算結果動態地調整了每一增量步上每一個積分點的黏聚力,從而實現了土體強度的各向異性,利用強度折減法求出抗滑樁加固各向異性邊坡的安全系數,并通過改變模型的樁間距,分析了土拱效應和安全系數在不同樁間距下的變化情況,得到了考慮了各向異性和水平土拱效應條件下的最優樁間距.

1 計算模型及基本算例

1.1 基本計算模型

本文基本計算模型取自文獻[23],為了充分發揮出抗滑樁的加固作用,令樁的埋深貫通土體,本文邊坡的寬度取2倍的樁間距以減小邊界效應的影響,具體參數如圖1所示.

圖1 基本模型計算尺寸

土體與抗滑樁的相關參數見表1.在建模的過程中,樁和土之間的切向接觸為Penalty,摩擦系數0.3,約為tan(0.75φ) ,法向接觸為Hard.整個模型采用底部固定、側面法向的約束形式,其他平面不設置約束條件.

表1 抗滑樁-邊坡體系強度參數

1.2 強度折減有限元法

在有限元軟件中,強度折減法根據迭代計算不收斂或者滑動面上的位移突變作為邊坡失穩依據,因此安全系數可視為使邊坡達到失穩狀態時的強度折減系數[24].通常根據下式進行計算：

式中：cm和φm為土體實際所需的抗剪強度；c和φ為土體所能夠提供的抗剪強度；Fs為強度折減系數.

1.3 基本算例及其結果驗證

基本算例中的土體設為均勻土體,未考慮各向異性,有限元模型網格如圖2所示,為了增加計算效率并保證計算結果的精確度,網格在滑動面與樁孔位置處進行了加密處理.

計算結果為：Fs=1.50,與原題1.536基本一致,相對誤差不足3%.微小的差別可能來自不同的網格劃分方法、網格的密度、增量步長等,通過對基本模型的驗證,說明了本文采用的建模過程是正確且可行的.

圖2 有限元網格

2 基于土體各向異性完善計算模型

2.1 各向異性在ABAQUS中的實現

在基本算例的基礎上,通過改變特定的土體參數,從而實現土體的各向異性.由于ABAQUS主程序只能將土體視為黏聚力與內摩擦角均為固定常數的各向同性土體,因此需要通過調用子程序實現土體強度和內摩擦角的變化.相對于黏聚力的各向異性,內摩擦角φ的各向異性產生的效果并不明顯,因此本文主要針對黏聚力的各向異性展開分析.各向異性的黏聚力通常按照下式確定[25]：

式中：cv和ch分別為垂直方向和水平方向的黏聚力；i為大主應力方向與垂直方向的夾角；k為各向異性系數,本文取k=0.5.

ABAQUS的編程通常是通過Fortran語言進行子程序的編寫而實現的,理想的各向異性需要通過UMAT子程序構建完整的摩爾-庫倫本構模型,本構模型的建立需要耗費巨大的工作量,因此本文根據王棟[22]采用的方法,使用主程序自帶的本構模型,通過USFLD子程序自定義不同的場變量對應的值(如圖3所示).

圖3 材料參數設置

將黏聚力按照大主應力方向分為若干小段,每一小段的黏聚力在該段內取均值,通過實時調用每個積分點上的應力分量計算大主應力方向,然后對大主應力方向大小的判斷確定場變量的序號,從而確定出該點的黏聚力取值,實現動態更新每個積分點的黏聚力.

2.2 考慮土體各向異性條件下的安全系數

通過強度折減法求解出安全系數,與原文結果對比發現,考慮了土體強度各向異性后,安全系數Fs為1.406,較未考慮土體各向異性時(Fs=1.536)有明顯的降低,未考慮土體各向異性時的計算結果會高估抗滑樁邊坡的穩定性.這是因為土體在考慮了各向異性條件時,黏聚力不再是某一固定的常數,而是根據i的不斷變化而改變,i的取值范圍為0°～90°,根據公式(5)可以發現,ci始終小于等于cv,因此當黏聚力整體呈現下降的趨勢時,土體的抗滑能力便會下降,滑坡的風險隨之上升,從而降低了安全系數.

3 樁間距對水平土拱效應以及安全系數的影響

3.1 樁間距對土拱范圍的影響

在土體各向異性的基礎上,對樁間距分別為2、4、8、12倍樁徑時的土拱范圍進行對比分析.由圖4可知,樁間土體的主應力方向由于樁內側與土體的摩擦而產生偏轉,在樁間和樁前也呈現出一定范圍的水平土拱效應,本文主要針對樁后的水平土拱效應研究樁間距的影響.如圖4所示,當樁間距逐漸增加時,樁后拱軸線的曲率逐漸降低,當s/d=12時,樁后水平土拱已經消失,說明隨著樁間距的增大,土拱越來越難以形成,當s/d=8以后,雖然能清晰地看到土拱的存在,但是其應力變化幅度已經非常小,可以理解為當s/d=8時,土拱效應已經非常不明顯.

圖4 水平土拱的變化與樁間距的關系

3.2 樁間距對中軸線上y方向的應力σy的影響

由圖5可知：當s/d≤6時,中軸線上的σy的突變較為明顯,x在-2～2 m的范圍內出現了兩個較為明顯的拱形,當s/d≥8時,中軸線上的σy為一條平緩的曲線,兩側沒有出現明顯的拱效應.因此,根據水平土拱效應的發揮程度可以確定出極限樁間距(使樁后能夠產生土拱效應的最大樁間距),由此可知模型的極限樁間距為s=8d.

圖5 中軸線上y方向應力分量與樁間距的關系

3.3 樁間距對y剖面上應力分量σx的影響

不同樁間距下的y剖面上應力分量σx的變化對比分析如圖6所示：離中軸線越近,應力分量越小,當樁間距在較小的范圍內(s/d≤3)時,水平土拱范圍也較小,并且樁前的水平土拱效應較樁后更加明顯；當樁間距逐漸增加時,樁后的土拱范圍逐漸增大至飽和的拱形；繼續增加樁間距,土拱范圍開始逐漸縮小.

圖6 樁前、后不同剖面上x方向應力分量

從應力變化曲線可以明顯地觀察到,當s/d=8時,拱形已經變形,從圓滑的曲線變成扁平的拋物線形狀,因此極限樁間距為8倍的樁徑.

3.4 樁間距對安全系數的影響

圖7 安全系數與樁間距的關系圖

圖8 最優樁間距的確定

圖7給出了不同樁間距下的安全系數曲線,如圖7所示：當樁間距逐漸增加時,安全系數呈現出逐漸降低的趨勢；當樁間距在較小的范圍內時,安全系數降低的速度較快,超過這一范圍時,繼續增加樁間距,安全系數降低的速度逐漸趨向于零.如果將曲線分為兩段,并將這兩段曲線各自擬合成一條虛直線,如圖8所示：虛線1的安全系數變化較大,而樁間距變化較小,因此虛線1以考慮安全系數為主；虛線2的安全系數變化較小,而樁間距變化較大,因此虛線2以考慮樁間距為主.兩條虛線的交點處既充分考慮了樁間距的影響,又充分保證了穩定性,故取兩條虛線交點的橫坐標為最優樁間距,約為s/d=5.

4 結 論

1)通過對ABAQUS的二次開發實現了土體的各向異性,證明了各向異性對抗滑樁加固邊坡有著不可忽視的影響,關于各向異性系數與初始黏聚力的大小對結果的影響還需進一步考究.

2)運用ABAQUS可以合理的模擬三維情況下土拱的形成過程,研究發現水平土拱的影響范圍不一定是隨著樁間距的增大而減小,在樁間距較小時,水平土拱的影響范圍隨著樁間距的增大而增大.

3)通過改變樁間距研究了不同路徑下的應力變化情況,分析了不同樁間距下安全系數的變化,較為準確地計算出極限樁間距以及最優樁間距的大小.

4)運用強度折減法和二分法可以較為準確地計算出土體在各向異性條件下的安全系數,但是手動輸入每一折減系數需要耗費較多的時間,因此在后續的研究中,安全系數的求解方法還需進一步優化.