凍融循環作用對土質邊坡穩定性的影響

宋彥琦， 鄭俊杰， 李向上， 李小龍， 邵志鑫

(中國礦業大學(北京)力學與建筑工程學院，北京 100083)

隨著季節的交替，凍土區會反復發生凍融循環作用，在此過程中，土壤理化性質的改變使邊坡穩定狀態逐漸遭受破壞，最終發生失穩。據調查，中國高海拔地區多年凍土面積居世界首位[1]，因此研究凍融循環作用對邊坡穩定性的影響具有重要且現實的意義。

中外對凍融循環作用下邊坡穩定性的研究取得了大量成果。實驗方面，常丹等[2]對青藏粉砂土進行了常規三軸實驗，結果表明隨著凍融循環次數的增加，土樣的彈性模量和靜強度先減小后增大，黏聚力逐漸減小，內摩擦角先減小后增加。楊更社等[3]對白堊系飽和砂巖進行了實驗研究，發現凍結作用對飽和砂巖的力學性質的影響較大。崔廣芹等[4]研究了內蒙古包頭凍土地區的粉砂土，發現對于不同干密度和不同含水率的土樣，隨著凍融次數的增加，黏聚力會減小，內摩擦角變化不大，黏聚力對抗剪強度影響較大。葛琪[5]對凍融界面的抗剪強度進行了研究，發現隨著凍融循環次數和含水量的增加，凍融界面土體的抗剪強度逐漸降低，且含水量的影響更大，但其影響隨凍融次數增加而遞減并最終穩定。梁波等[6]對青藏鐵路清水河、北麓河的地基土進行取樣進行室內模擬試驗，發現隨著凍融循環次數的增加，土體密度逐漸趨于某一定值，融沉總變形在增加的同時呈現指數遞減的趨勢。通過這些實驗可以看出，凍融循環作用是影響巖土材料強度的一個重要因素。數值模擬方面，鄭穎人等[7-8]和趙尚毅等[9]通過有限元軟件利用強度折減法對邊坡安全系數進行求解，發現有限元強度折減法完全可以用來解決工程實際問題，并對有限元強度折減法的計算精度和影響因素進行了探討和完善。張永明等[10]和陳印東等[11]通過強度折減法利用有限元軟件對邊坡進行模擬，其計算結果與極限平衡法基本一致，表明強度折減法在邊坡穩定性的判定上具有合理性與可行性。史俊濤等[12]利用容重增加法和強度折減法對某土質邊坡進行穩定性分析，發現強度折減法相比容重增加法使用更加合理。

綜合上述研究可以看出，雖然中外在凍融循環作用對巖土力學性能的實驗研究以及邊坡穩定性的數值模擬研究已經取得了一定的成果，但針對凍融循環這個影響因素對邊坡整體穩定性的數值模擬研究很少，且強度折減法對于邊坡穩定性的評價具有合理性和實用性。本文利用強度折減法，通過大型有限元軟件ABAQUS計算獲得邊坡的穩定安全系數，與簡化Bishop法獲得的理論安全系數進行比較，并以特征部位位移發生突變作為破壞判據，得出了邊坡安全系數隨坡角的減小而增大，隨凍融循環次數增加而逐漸降低最終趨于平穩的變化趨勢，該文可為邊坡穩定性影響因素的探究提供參考。

1 邊坡穩定性的數值模擬

1.1 強度折減法的基本原理

強度折減法是在理想彈塑性的有限元計算中，逐漸降低邊坡巖土體的抗剪強度參數，直至達到極限破壞狀態，同時得到邊坡的強度儲備安全系數Fs[13]。

對于摩爾-庫侖材料，關于剪應力τ有：

(1)

所以有：

(2)

(3)

式中，σ為法向應力，c和φ為有效黏聚力和有效內摩擦角；c′和φ′為折減后的有效黏聚力和有效內摩擦角。

目前有限元中邊坡整體失穩的判據主要有三種[8]：①以有限元計算是否收斂作為破壞判據；②以特征部位位移發生突變作為破壞判據；③以是否形成連續的貫通區作為破壞判據。

1.2 計算模型和材料參數

某邊坡模型尺寸如圖1所示，土樣的力學參數采用文獻[14]中尾礦土的數據，其中土體密度ρ=1.492 g·cm-3，泊松比取v=0.35，彈性模量E采取圍壓σ3=0.2 MPa時測得的數值，不同凍融循環次數下土體物理力學參數如表1所示。

圖1 模型尺寸Fig.1 Model size

表1 不同凍融循環次數下土體物理力學參數

按照二維平面應變建立有限元模型，計算采用Mohr-Coulomb準則，模型邊界條件為：左右兩側限制水平方向的位移，底部限制水平和垂直兩個方向的位移，上部為自由邊界。模型只受重力荷載作用，整個模型包括349個單元，單位類型為CPE4，網格劃分如圖2所示。

圖2 有限元網格劃分Fig.2 Finite element mesh

1.3 計算結果

1.3.1 塑性應變分析

塑性應變等值云圖反映了邊坡塑性應變的分布，圖3所示為1次凍融循環作用下不同折減系數時三種不同情況下邊坡的塑性應變等值云圖，邊坡變形顯示比例設置為0。

圖3 不同折減系數下的塑性區Fig.3 Continuous contour of plastic strain under different reduction factors

由圖3塑性應變的分布可直觀判斷出邊坡失穩的發展過程，即塑性區最先從坡腳處產生，隨著邊坡發生屈服的區域逐漸擴大并向坡頂延伸，直至到達坡頂，即坡腳與坡頂間的塑性區發生貫通，邊坡發生失穩。

1.3.2 滑動面確定

滑動面的位置與形狀可以通過位移等值線分布圖來確定。

圖4為1次凍融循環作用下有限元計算不收斂時的位移等值云圖。圖4清楚地定位了滑動面的位置，形狀呈大致的圓弧狀。

圖4 用位移等值云圖表示的滑動面位置和形狀Fig.4 Position and shape of the sliding surface represented by the displacement equivalent cloud

1.3.3 安全系數分析

在計算過程中，將場變量取為強度折減系數，在0.5～7變化，線性增加場變量的大小，增量為0.25，特征點取為邊坡坡面左上角頂點，見圖1紅點。圖5為1次凍融循環作用下特征點處折減系數FV隨水平位移U1的變化關系圖。

圖5 折減系數隨位移的變化關系Fig.5 Variation of the reduction factor with displacement

如圖5所示，A(-0.001 3,4.954 2)處特征點位移發生突變，B(-0.629 2,5.044 3)處計算不收斂。折減系數從0.5～4.954 2這個階段，特征點位移幾乎一直為0，說明這個階段坡頂面幾乎沒有產生任何變形，由圖3(a)～圖3(c)可看出當折減系數未達到4.954 2時塑性區并未蔓延至特征點，因此在這個階段，邊坡尚處于穩定階段。直至折減系數大于4.954 2之后，水平位移絕對值迅速增加，說明這個階段特征點周圍土體受到了很大的影響，導致坡頂面迅速產生大量變形，位移為負值說明該部位處于向下滑動的變形狀態，由圖3(d)可看出當折減系數達到4.954 2時塑性區延伸至特征點，邊坡產生了一個完整的塑性貫通區，當位移達到一定值時，邊坡發生失穩。

1.3.4 安全系數確定

建立相同的邊坡計算模型，如圖1和圖2所示，分別采用凍融循環次數為0、1、5、10、15次時土體的力學參數，如表1所示，計算邊坡的安全系數，選取不同判據得出的安全系數如表2所示。其中，安全系數F1、F2、F3分別為以特征部位位移發生突變作為破壞判據、以是否形成連續的貫通區作為破壞判據、以有限元計算是否收斂作為破壞判據所得出的結果。

表2 不同凍融循環次數下不同破壞判據求出的邊坡安全系數

由表2可看出采用有限元中邊坡整體失穩的三種判據得到的安全系數相差不大，其中以特征部位位移發生突變作為破壞判據與以是否形成連續的貫通區為破壞判據得出的安全系數更為接近。

2 邊坡穩定性的理論分析

2.1 邊坡穩定分析的簡化Bishop法

采用文獻[15]中簡化Bishop的積分表達求解邊坡的安全系數。簡化Bishop的積分表達具體如下：邊坡高度為H，坡角為β，邊坡危險滑移線為圓弧，以坡角為坐標原點，水平方向為x軸，垂直方向為y軸，滑移線的圓心為O(a,b)(a，b為圓心在x，y軸的坐標)，半徑為R，將邊坡滑體部分分為若干條塊，條塊寬度為積分變量dx，條塊高度h的表達式為

(4)

式(4)中：

邊坡安全系數F為

(5)

式(5)中：γ為滑體的容重，α為內摩擦角，c為有效黏聚力。

2.2 滑動圓弧圓心范圍的確定

簡單土質邊坡最危險圓弧圓心的范圍，可按潘家錚[16]的方法確定，該方法如圖6所示。

圖6 滑動圓弧圓心范圍的確定Fig.6 Determination of the center of the sliding arc center

邊坡坡度為1：m，高為H，長為L=mH，作出邊坡的中垂線和中法線，以邊坡中點O為圓心，以R1=L/2和R2=3L/4為半徑作圓弧，在邊坡上側分別與中垂線和中法線圍成一個四邊形，最危險圓心常常位于這個四邊形范圍內。

在搜索最危險圓心時，為方便編程，將圓心范圍擴大為四邊形的外接矩形CDEF，如圖6所示。

2.3 計算結果

通過MATLAB編程實現邊坡的安全系數的求解，其中計算步長設置為0.1 m，計算結果如表3所示。

表3 不同凍融循環次數下邊坡安全系數的理論解

3 結果分析

3.1 不同破壞判據與理論解的對比分析

表4為凍融循環次數分別為0、1、5、10、15次時，邊坡采取強度折減法選取不同破壞判據和簡化Bishop法得出的理論解的誤差對比，不同破壞判據結果的誤差取5次凍融循環情況的誤差平均值。其中破壞判據一為特征部位位移發生突變，破壞判據

表4 不同破壞判據求出的安全系數與理論解的誤差對比

二為形成連續的貫通區，破壞判據三為有限元計算不收斂。

由表4可看出，分別采用三種破壞判據得出的不同凍融循環條件下邊坡的安全系數與簡化Bishop法得到的理論值的誤差都很小，即通過有限元模擬與簡化Bishop法得到的安全系數基本一致，說明這兩種方法在本例中都是可行的。

通過判據一得出的安全系數與理論解的誤差和通過判據二得出的安全系數與理論解的誤差相差甚小，且都小于通過判據三得出的安全系數與理論解的誤差，即以特征部位位移發生突變和是否形成連續的貫通區作為破壞判據相比以有限元是否收斂作為破壞判據更為準確，分析其原因是當塑形變形區從坡腳貫通到坡頂時，若無其他約束條件比如抗滑樁、錨桿等邊坡支護設施，塑性區上方部位會迅速產生大量變形，位移自然會大幅度增加，但此時由于載荷步步長設置或其他原因導致有限元計算不一定不收斂，即有限元計算收斂時不一定表明邊坡處于穩定狀態。因此在實際工程應用中，可根據實際情況，將特征部位位移是否發生突變和塑性區是否貫通兩個判據綜合考慮判斷邊坡的穩定性狀態。

在接下來對邊坡安全系數的求解過程中，通過有限元強度折減法進行計算且全部以特征部位位移發生突變為破壞判據確定安全系數數值。

3.2 凍融循環次數對邊坡穩定性的影響

為更直觀地分析得到凍融循環次數對邊坡安全系數的影響，將表2中以特征部位位移發生突變為破壞判據得到的安全系數繪成折線圖，如圖7所示。

圖7 安全系數隨凍融循環次數的變化關系Fig.7 Variation of safety factor with the number of freeze-thaw cycles

由圖7可看出，在本例中，當坡比固定時，凍融循環作用后，邊坡穩定性逐漸下降，從10次凍融循環到15次凍融循環，有限元計算的邊坡的安全系數下降比僅為7.494 3%，而經1次凍融循環作用后，安全系數下降比達24.959 1%，下降幅度最大。

3.3 邊坡坡比對邊坡穩定性的影響

為研究邊坡坡比對邊坡穩定性的影響，分別選取坡比為1∶0.5、1∶0.75、1∶1、1∶1.25建立不同坡比的邊坡模型，如圖8所示，計算相應邊坡在凍融循環次數分別為0、1、5、10、15次時的穩定系數，結果如圖9所示。

圖8 不同坡比邊坡尺寸示意Fig.8 Slope size at different slope ratios

圖9 安全系數隨坡比的變化關系Fig.9 Variation of safety factor with slope ratio

由圖9可清楚地看出，當坡比一定時，隨著凍融循環次數的增加，四種不同坡比情況下的邊坡的安全系數下降的幅度都會逐漸減小，并最終趨于穩定；當凍融次數一定時，四種不同坡比情況下的邊坡的安全系數都會隨著坡角的減小逐漸增大，即邊坡越平緩，安全系數越大，邊坡穩定性越好。

4 結論

(1)通過強度折減法，采用有限元中邊坡整體失穩的三種判據求出一定坡比下不同凍融循環次數下的安全系數，三種判據得到的結果相差不大。將其與簡化Bishop法得到的理論安全系數進行對比，兩者計算結果基本一致。

(2)以特征部位位移是否發生突變和是否形成連續的貫通區作為破壞判據相比以有限元是否收斂作為破壞判據更為準確，在實際工程應用中，可根據實際情況，將特征部位位移是否發生突變和塑性區是否貫通兩個判據綜合考慮判斷邊坡的穩定性狀態。

(3)通過強度折減有限元方法建立ABAQUS力學模型，研究邊坡的穩定性，相比簡化Bishop法而言，更加客觀簡便，不需要假定危險滑動面，就能夠直接動態地展示不同邊坡的失穩破壞過程，該方法可為邊坡穩定性影響因素的探究提供參考。

(4)討論了凍融循環以及邊坡波比對邊坡穩定性的影響，但是邊坡研究還應考慮現場的環境溫度，濕度等因素的影響，另外力學模型的簡化程度等因素對于模擬結果都是有影響的，這些都是今后研究中需要進一步考慮的問題。