各向異性成層邊坡的穩定性分析

羅凌暉，周建，蔡露，溫曉貴

(1.浙江大學濱海和城市巖土工程研究中心，浙江杭州，310058；2.浙江省城市地下空間開發工程技術研究中心，浙江杭州，310058；3.蘇州常宏建筑設計研究院姑蘇分公司，江蘇蘇州，215000)

邊坡穩定分析是邊坡設計的前提，能有效判斷邊坡失穩時的推力，為后續支護結構的設計提供依據，從而有效防止因邊坡失穩破壞而形成的各種地質災害(如滑坡、崩塌、泥石流等)，為構筑物的正常使用和人們的生命財產安全保駕護航[1]。常用的邊坡穩定分析方法有極限平衡法、極限分析法、滑移線場法、有限元數值分析方法等等。人們分析時往往將邊坡土體視為各向同性材料[2-5]，而天然軟土由于沉積，初始不等向固結及復雜應力路徑等原因呈現出顯著的各向異性。ZAREI等[6-11]的研究結果表明：主固結方向土體不排水抗剪強度要比其他方向的高。對于開挖邊坡而言，邊坡圓弧滑動面上土體單元的大主應力軸方向角是連續變化的，滑動面底部區域土體單元的大主應力軸與豎直方向夾角α甚至可以達到90°，由于各向異性的客觀存在，該區域土體強度低于邊坡頂部土體強度，如果統一采用三軸壓縮不排水抗剪強度指標對土坡進行穩定分析，可能會得到偏高的安全系數。因此，如果不考慮土體各向異性，那么在對巖土工程問題的分析中將無法保證結果的準確性[12-13]。目前，考慮土體強度各向異性的邊坡穩定的主要研究方法有傳統分析方法和理論分析公式法。傳統分析方法首先假設有1 個確定的滑動面。例如，CHEN 等[14]將Casagrand 土體強度各向異性表達式引入極限分析法的上限理論，得出了不同坡度角以及不同各向異性系數條件下的邊坡穩定系數；趙紅華[15]采用Bishop條分法分析了土體強度各向異性對邊坡安全系數的影響，得出土體黏聚力各向異性對邊坡穩定性的影響要大于內摩擦角各向異性的影響，且各向異性程度愈大，影響愈強烈等。理論分析公式法僅僅是三軸壓縮抗剪強度與三軸拉伸抗剪強度的簡單組合，沒有考慮到土體實際的力學特性。例如，王棟等[16]以Casagrand 公式為基礎開發了能自動搜索邊坡安全系數的有限元數值分析方法，并探究了邊坡安全系數隨各向異性系數的變化規律；SU 等[17]重新建立了以三軸壓縮和三軸拉伸抗剪強度為變量的土體強度各向異性表達式，分析了強度各向異性對邊坡安全系數的影響。然而，將試驗所得土體強度各向異性引入商業有限元軟件中從而對各向異性成層邊坡進行穩定分析的研究較少。本文作者以杭州軟黏土復雜應力路徑下不排水抗剪強度試驗及屈服特性試驗結果為基礎，提出考慮土體不排水抗剪強度各向異性的半經驗公式，并依據該公式的基本原理提出考慮土體強度各向異性的修正Von-Mises屈服準則；利用商業有限元軟件ABAQUS 的UMAT子程序進行二次開發，建立考慮土體強度各向異性成層邊坡分析模型；在此基礎上，采用強度折減法探究成層邊坡土體各向異性和不排水抗剪強度等因素對邊坡安全系數的影響，并結合實際算例，驗證本文方法的可靠性。

1 各向異性不排水抗剪強度公式推導

關于不排水抗剪強度的研究，大多是通過與沉積面呈不同角度切取土體制作試樣，然后，利用這些試樣進行常規室內土工試驗，如常規三軸、直剪和平面應變試驗等。此類試驗由于儀器的局限性，無法真正實現主應力軸的旋轉[12-13]。為了真正實現土樣加載直至破壞過程中主應力軸方向的偏轉，沈揚[9]應用空心圓柱扭剪儀(HCA)探究了杭州地區軟黏土強度各向異性，發現其不排水抗剪強度各向異性隨主應力軸方向呈現“勺型變化”。高彥斌等[10]利用室內異形十字板試驗對上海軟黏土強度各向異性特征進行了研究，發現試驗結果與Casagrande 公式和Carrillo 公式[18]計算結果較吻合。應宏偉等[19]結合Casagrande公式中考慮不排水抗剪強度最小和最大時所對應的α，對文獻[9]中的試驗結果進行了擬合，得到杭州地區典型軟黏土在大主應力方向為α剪切時的不排水抗剪強度Suα為

式中：Suv為三軸壓縮不排水抗剪強度；α為大主應力軸方向與豎直方向夾角，取值范圍為0°～90°；k為各向異性比；根據空心圓柱扭剪儀(HCA)試驗中空心土樣的雙向軸對稱特性，α為-90°～0°時對應的Suα可由0°～90°時的Suα對稱得到；k取值為α= 90°和α= 0°時不排水抗剪強度的比值。

主應力軸方向常常和破壞面的方向聯系在一起，LO[20]的研究結果表明，主應力軸方向與土體剪切破壞面的夾角ψ為常量。

式中：?′為土體的有效內摩擦角。沈揚[9]根據主應力軸方向旋轉的不排水抗剪強度試驗，發現剪切破壞面與主應力軸方向的夾角隨主應力軸方向改變不明顯，一般為25°～30°。杭州地區典型軟土為全新世最后一次海侵過程中形成的海相沉積土，呈現流塑狀態并具有水平層理。王冠英等[21]通過微觀電子顯微鏡對杭州地區典型軟土進行定向度測定，發現該土樣顆粒薄平表面與沉積面平行的顆粒數約為與沉積面垂直的顆粒數的2 倍。由此可見，當剪切破壞面方向為水平向時，土體顆粒之間最易產生相對滑動。

本文根據上述研究成果以及杭州地區典型軟土的微觀破壞機理，在α=?′/2 - 45°方向上取不排水抗剪強度為最大值，在α=?′/2 + 45°方向上取為最小值，即土體剪切破壞面為豎直方向和水平方向。根據Casagrande 公式，則有

式中：Suh為三軸拉伸不排水抗剪強度。

將式(3)中的大主應力軸方向角轉換為剪切方向角可以得到

將α= 0°和α= 90°分別代入式(4)整理后可得：

由式(5)和式(6)可得：

由式(7)和式(8)可得：

最后，將式(9)和式(10)代入式(4)得到不排水抗剪強度表達式：

式(11)即為考慮土體有效內摩擦角影響的強度各向異性半經驗公式，包含三軸壓縮不排水抗剪強度Suv、三軸拉伸不排水抗剪強度Suh以及土體有效內摩擦角?′共3個參數。

為了探究土體各向異性程度對巖土工程問題各種系數的影響，通常將Suh/Suv定義為一個新的變量，但尚未形成一致的表達[16-19]。本文統一將k=Suh/Suv定義為不排水抗剪強度各向異性比，簡稱為各向異性比。將k代入式(11)可以得到

對比式(1)和式(12)可知：式(1)是式(12)取?′=30°時的一種特殊情況，因此，本文公式適用性更廣。圖1所示為在不同情況下，大主應力軸方向角-90°≤α≤90°時，歸一化不排水抗剪強度Suα/Suv隨α的變化情況。圖1中，本文公式取有效內摩擦角?′=30°。

由圖1可知：本文公式與Casagrande公式的強度各向異性規律具有顯著差異，本文公式可以考慮大主應力軸方向角在-90°≤α≤90°范圍內的變化規律，最大值并不在α= 0°處，而是在-90°≤α≤0°范圍內，這與文獻[10]中的試驗結果相符。相比于Casagrande 公式，本文公式歸一化抗剪強度Suα/Suv在0°≤α≤90°范圍內的結果與文獻[9-10]中的試驗結果擬合度更高，可見本文公式在參數選取合理的前提下具有一定的優勢。因此，下面將以式(12)為標準，推導考慮強度各向異性的修正Von-Mises準則。

圖1 不同公式及試驗所得歸一化不排水抗剪強度Suα/Suv對比Fig.1 Comparison of normalized undrained shear strength Suα/Suv by different formulas and various test results

2 考慮強度各向異性的修正Von-Mises屈服準則

在工程實踐中，對基坑邊坡開挖進行不排水數值分析時，通常選用固結不排水抗剪強度指標ccu和?cu，按照Mohr-Coulomb 或Drucker-Prager 屈服準則進行分析。對于海底邊坡、水庫邊坡或地下水位較高的基坑開挖邊坡，土體主要由軟黏土構成且處于飽和狀態，滲透系數較小，加荷速率較快。龔曉南等[22-23]認為應選用不排水抗剪強度指標對該類問題進行計算分析。不排水條件下的軟黏土力學特性更符合Tresca或Von-Mises 屈服準則。由于在三維主應力空間中Tresca 屈服面有角點存在，在數值分析中會遇到困難，因此，本文采用Von-Mises屈服準則和強度各向異性半經驗公式推導考慮強度各向異性的修正Von-Mises屈服準則。

為了計算方便，巖土工程中常用有效主應力組成的應力不變量來表示土體的應力狀態。常用的應力不變量包括平均有效應力p′、等效剪應力J和應力羅德角θ。

式中：σ′1，σ′2和σ′3分別為大主應力方向、中主應力方向和小主應力方向的有效正應力。則Tresca準則的屈服函數fT可以用主應力表示為

式中：Su為土體不排水抗剪強度。

Von-Mises 準則的屈服函數fM可以用主應力表示為

式中：q為等效應力，kf為Von-Mises 準則的屈服強度。將式(16)和(17)用應力不變量表示為：

Tresca 準則和Von-Mises 準則的屈服函數在偏應力平面內的形狀分別為正六邊形和圓形，有2個Von-Mises 圓，其中，一個是Tresca 正六邊形的外接圓，另一個是內切圓。一般認為Von-Mises 屈服函數是Tresca 屈服函數的近似形式，因此，需要尋找和Tresca 正六邊形擬合最好的圓。由文獻[24]可知擬合最好的圓在θ=±15°處，因此，將θ=±15°代入式(18)，聯立式(18)和式(19)并取fT=fM，可以得到

為進一步考慮不排水抗剪強度各向異性的影響，將式(20)中的Su和式(12)中的Suα對應起來，將式(12)和式(20)代入式(17)，得到修正的Von-Mises 準則的屈服函數為

修正的Von-Mises準則仍假定塑性應變增量與塑性勢面正交，則土體不會產生塑性體積應變增量，以此來保證模型能準確反映飽和土體的不排水力學特性。

為了完整表達該模型，還需要確定彈性參數。由于不排水剪切過程中不會產生彈性體積應變，取泊松比υ≈0.5，這樣，模型中就只剩不排水彈性模量Eu和不排水剪切強度Su需要提前確定。在實際工程應用中，由于土并非理想彈塑性體，其變形包含了可恢復的彈性變形和不可恢復的殘余變形2 個部分，因此，在靜荷載作用下計算土體變形時常采用壓縮模量和變形模量等。本文采用修正的Von-Mises 屈服準則，將土體視為彈性-理想塑性材料，因此，對于不排水彈性模量這一參數應選用土體彈性模量進行計算分析。一般通過室內三軸儀進行不排水三軸壓縮試驗，根據得到的應力-應變關系曲線確定土的彈性模量。詳細的試驗方法參見文獻[25]。

3 強度折減法有限元模擬

隨著計算機應用技術的高速發展，利用有限元強度折減法進行邊坡穩定分析受到國內外研究者的廣泛關注[26-28]。相較于傳統的邊坡穩定分析方法如極限平衡法、極限分析法和滑移線場法等，有限元強度折減法不需要事先假定邊坡最危險滑動面的形狀和位置，只需要不斷降低邊坡土體的強度參數，進而使邊坡土體因抗剪強度不足而發生失穩破壞，就能得到邊坡的滑動面形狀及相應的安全系數。

目前，對于均質邊坡，有限元強度折減法已經取得了較可信的結果。該方法的核心理念如下：在理想彈塑性有限元計算中逐漸降低邊坡巖土體抗剪切強度，直到其達到破壞狀態為止[4]，此時，得到的強度折減系數即為邊坡穩定安全系數F，相關原理可以表示為

式中：c和φ為土體原本的抗剪強度參數；cm和φm為折減后土體實際用于計算的抗剪強度。

根據有限元分析結果判定邊坡處于破壞狀態，是運用強度折減法求解安全系數的重點。根據目前的研究進展，邊坡失穩判據可分為如下3種：1)有限元程序迭代過程中出現力或位移的不收斂，即最終數值計算結果不收斂；2)特征點發生位移突變或具有較大的變形趨勢；3)從坡腳到坡頂形成等效塑性應變或廣義塑性應變的連續貫通區。

根據文獻[28]，將邊坡滑動面上所有單元節點中任一節點的塑性應變或位移出現突變且產生較大程度無限制的塑性流動作為邊坡失穩的標志。每一折減系數最多允許500次迭代，采用平面應變四邊形八節點等參單元(CPE8)對各向異性成層邊坡模型進行網格劃分。

為了模擬土體強度各向異性，利用ABAQUS 的UMAT子程序編寫修正的Von-Mises屈服準則，計算流程如下。

1)體力施加采用分級加載，每一級加載開始時，主程序調用UMAT 子程序的材料信息傳入單元積分點，同時傳入時間步長、應變增量、荷載增量、當前狀態應力和應變以及其他所有求解過程所需要的變量。

2)UMAT 子程序沿著應變增量路徑對本構方程進行積分，得到新的應力增量，并同時更新其他相關變量，然后，提供Jacobian矩陣給主程序以計算形成新的整體剛度矩陣。

3)ABAQUS 主程序結合當前加載步的殘余荷載和位移增量進行平衡檢驗，不平衡力和位移增量修正容許的相對誤差分別為0.5%和1.0%，若不滿足平衡收斂條件，則繼續迭代，直至收斂[16]。然后，進行下一加載步計算。

本文采用完全隱式的向后Euler 積分回映算法進行計算。該算法的核心是將增量彈塑性本構方程轉化為一組非線性方程組，然后，利用牛頓迭代方法對該非線性方程組進行線性化處理并求解。隱式回退算法因強化了子增量步結束時的一致性，避免了屈服面的漂移，具有較高的精確性。

4 算例驗證與討論

由于很難假定統一的滑動面形式，傳統邊坡穩定分析方法難以推廣到各向異性成層邊坡，而有限元方法在計算完成后自動得到滑動面，適用于確定實際工程中廣泛存在的復雜成層邊坡的安全系數。下面根據文獻[29]中的算例，驗證本文方法對各向異性成層邊坡進行有限元計算分析的適用性。

4.1 算例驗證

成層邊坡示意圖如圖2所示。成層邊坡均由不排水摩擦角φu= 0的黏土構成，土體容重γ= 20 kN/m3，有效內摩擦角φ′= 20°，上層邊坡土體的豎直向黏聚力Su1=25 kPa，將下層邊坡土體的豎直向黏聚力Su2與上層豎直向黏聚力Su1之比定義為強度比。應用不排水指標分析，取泊松比υ≈0.5，Eu= 400Su(與文獻[29]中算例的彈性模量相等)。采用本文第2節中提出的相關聯的修正Von-Mises準則，利用強度折減法進行有限元分析計算。邊界條件設置為邊坡上邊界自由約束，左右邊界水平約束，底部完全約束。

圖2 成層邊坡示意圖Fig.2 Diagram of layered slope

GRIFFITHS 等[29]在計算時沒有考慮土體強度各向異性，因此，這里取各向異性比k= 1 進行分析。經ABAQUS有限元軟件計算后，結合文獻[29]中的計算結果，得到邊坡安全系數與強度比Su2/Su1的關系如圖3所示。由圖3可以看出：本文擬合結果與文獻[29]中的計算結果具有相同的變化趨勢，本文采用的是修正Von-Mises 準則，因此，本文結果較采用Mohr-Coulomb 屈服準則計算所得邊坡安全系數大3%～7%。

圖3 邊坡安全系數與強度比(Su2/Su1)的關系Fig.3 Relationship between safety factor and strength ratio(Su2/Su1)

由圖3還可以看出：邊坡安全系數與強度比曲線存在1個明顯的轉折點，即強度比Su2/Su1約為1.5的位置。根據Su2/Su1的不同，提取ABAQUS 分析計算后3種典型的邊坡失穩破壞模式，如圖4所示。

分析圖3和圖4可知：當Su2/Su1＜1.5時，邊坡失穩破壞的滑動面呈現典型的坡底圓形式，在這一階段，隨著強度比Su2/Su1逐漸增大，邊坡安全系數基本呈現出線性增長模式；當Su2/Su1＞1.5 時，邊坡失穩破壞的滑動面呈現典型的坡腳圓形式，此時，由于邊坡上、下層土體的抗剪強度差異較大，當上層土體因強度折減而發生失穩破壞時，下層土仍具有較大的抗剪強度，還未產生較大的塑性應變，因此，這一階段所對應的安全系數幾乎保持不變；當Su2/Su1≈1.5，有限元計算不收斂時，2種失穩破壞形式同時體現在這個邊坡上(見圖4(b))。

圖4 不同強度比對應的邊坡破壞形式Fig.4 Different failure types of slope corresponding to different Su2/Su1

同時，根據圖4中的邊坡破壞形式可以發現：隨著Su2/Su1不斷增大，邊坡的失穩破壞形式逐步由坡底圓向坡腳圓發展，且當Su2/Su1約為1.5 時，邊坡顯示出2 種破壞形式都有所體現的疊合狀態，這與文獻[29]中的有限元分析結果一致。不同文獻中成層邊坡穩定性分析結果與本文計算結果對比如表1所示。由于不同文獻中算例所采用土體參數均不一致，因此，沒有選擇對土體參數較為敏感的安全系數作為比較對象，僅從破壞形式方面與本文計算結果進行對比。由表1可以看出本文計算結果與文獻[16，30-34]中的計算結果一致。

綜上可知，本文的計算結果較為合理。下面針對以上3種典型的成層邊坡破壞模式探究土體強度各向異性對安全系數的影響。

4.2 土體強度各向異性的影響程度

成層邊坡算例的安全系數隨各向異性比和強度比的變化如表2所示。由表2可知：不同強度比下邊坡安全系數隨各向異性比k的變化呈現相同的趨勢。值得注意的是，當取各向異性比k= 0.8 時，安全系數反而較各向同性邊坡(k= 1.0)的有所提升。這是因為，當利用本文式(12)進行分析計算時，不排水抗剪強度最大值和最小值并不在α= 0°和α= 90°方向上，而是出現在α=?′/2 - 45°和α=?′/2 + 45°方向上。因此，當各向異性比接近于1.0 時，若根據峰值抗剪強度確定成層邊坡安全系數，則可能出現考慮強度各向異性后，邊坡安全系數反而提高的現象。

為了直觀地分析各向異性比對安全系數的影響，以各向同性成層邊坡計算得到的安全系數為基準，得到邊坡安全系數變化幅度T與各向異性比k的關系，如圖5所示。

結合表2和圖5可知：整體而言，當各向異性比較小(k≤0.8)時，無論Su2/Su1如何取值，隨著各向異性比不斷增大，邊坡安全系數將逐漸增大，且當土體呈現較強的各向異性時(k= 0.4)，在3 種強度比下，有限元計算所得安全系數比各向同性邊坡安全系數最大可減少50.77%；但對于不同強度比的成層邊坡，其影響程度卻呈現出一定的差異性。隨著各向異性比增大，當Su2/Su1=1.0 即對應滑動面呈現典型的坡底圓形式時，安全系數首先迅速增大，之后增大趨勢逐漸平緩；當Su2/Su1=2.0 即邊坡最危險滑動面為坡腳圓時，安全系數呈現出先平緩增大后迅猛增大的趨勢；而當Su2/Su1=1.5時，安全系數近似呈線性增大。這是因為：當各向異性比較小時，Su2/Su1=1.5 的最危險滑動面更傾向于坡底圓形式；隨著各向異性比增大，其最危險滑動面向坡腳圓形式發展，直至k= 1.0 時對應的成層邊坡失穩破壞呈現出2種滑動面形式都有所體現的疊合態。

實際工程中的邊坡土體通常都存在不同程度的強度各向異性，但目前利用傳統邊坡穩定分析方法或有限元法評價其穩定性時大多按各向同性邊坡處理，尤其是對于成層邊坡，若不考慮各向異性比的影響，則往往容易過高估計成層邊坡的穩定性，各向異性程度越大，相對誤差越大。

表1 成層邊坡穩定分析計算結果對比Table1 Comparison of computing results of layered slope stability analysis

表2 安全系數隨各向異性比和強度比的變化Table2 Variation of safety factor with the anisotropy ratio and strength ratio

圖5 安全系數變化幅度T與各向異性比k的關系Fig.5 Relationship between safety factor change amplitude and anisotropy ratio

5 結論

1)提出能考慮土體有效內摩擦角的影響的各向異性不排水抗剪強度公式，同時得到能考慮土體強度各向異性的修正Von-Mises屈服準則，并結合有限元強度折減法編寫了能計算各向異性成層邊坡安全系數的ABAQUS 有限元模型；相較于傳統邊坡穩定分析方法，本文方法更為簡便、準確，且當土體參數選用得當時，即使選用簡單的彈塑性模型進行計算分析，同樣可以得到較合理的計算結果。

2)當強度比Su2/Su1不同時，成層邊坡的最危險滑動面呈現出3種不同的破壞模式，無論是何種破壞模式，若不考慮土體各向異性影響，則安全系數計算結果不準確；當土體呈現較強的各向異性時(k= 0.4)，有限元計算得到的安全系數甚至能比各向同性邊坡的安全系數減少約50%。