多年凍巖土區露天礦邊坡局部穩定性探究

李國鋒，李 寧,2，劉乃飛，朱才輝

(1. 西安理工大學 巖土工程研究所，陜西西安710048；2. 蘭州交通大學土木工程學院，甘肅蘭州730070； 3.西安建筑科技大學土木工程學院，陜西西安710055)

多年凍土和季節凍土區的面積約占地球陸地面積的50%，我國永久性凍土和季節性凍土的分布區域占國土面積的70%以上[1]。青海木里煤礦就處在青藏高原高寒帶山地的多年凍土區，是青海省最大的煤炭資源分區之一。多年凍土區的露天礦邊坡有其特殊的地層分布，本文對王曉東等[2]提出的地層概化模型進行優化，如圖1所示。

圖1 多年凍土區露天礦邊坡的地層特性Fig.1 Strata status of open-pit slope in permafrost regions

永凍層上部為受溫度影響的季凍層，邊坡開挖后，坡表一定深度內的永凍層會變成季凍層，且永凍層下的不凍層在坡表處也會出現凍融交替現象。

該地區年最低氣溫-34 ℃，最高氣溫19.8 ℃，年平均氣溫為-4.2 ℃～-5.1 ℃，年平均地溫-1.0 ℃～-3.5 ℃。由于氣候的季節性冷暖變化，季凍層每年4月下旬開始融化，9月底至10月中旬達到最大融深，融化速度各月稍有差異，平均0.01 m/d，9月末氣溫開始下降，當進入負溫后，季融層自上而下迅速回凍，至12月初季融層全部凍結。即在氣溫下降回凍開始或回凍后一段時間內會出現最大融化深度。張寶龍[3]、張宏剛等[4]經近 2 年多觀測，認為該處季凍層厚約5～6 m，多年凍土層厚約58～60 m。徐拴海[5]進一步指出，季凍層在表層5 m 左右波動，并基于調查研究，首次提出了凍融巖土邊坡5種基本破壞模式，指出逆層邊坡多為中上部局部突出而下部崩塌，順層邊坡多為層間塊片狀脫落滑塌，切層邊坡多為巖石劣化塊體滑落或大塊傾倒，破碎地層主要為表層剝蝕，而松散地層為掏蝕、滑塌等，并明確指出環境的反復劇烈變化主要影響坡體淺層地表，而其中溫度變化是巖體性態改變的驅動力，水分是巖體凍融劣化的根源。

目前，關于多年凍土區的露天礦邊坡穩定性分析評價方法，主要集中于強度折減法、極限平衡法及有效應力法。武鶴等[6]采用極限平衡法推導了寒區土質邊坡凍融滑塌安全系數表達式，并討論了各項參數對邊坡穩定性的影響。曹蘭柱等[7]，基于強度折減法探討了不同開挖位置、開挖坡角和內排追蹤距離下的邊坡變形破壞機制和其對穩定性的影響規律。李榮建等[8]研究表明，傳遞系數法在分析邊坡的局部滑動問題中具有明顯的局限性。侯玲等[9]基于場變量的邊坡有限元強度折減法實現了高效邊坡安全系數求解，并基于最大等效塑性應變原理實現了邊坡潛在滑動面位置的精細搜索。劉紅軍等[10]基于有效應力法推導了季節性凍土區邊坡穩定性分析方法，且認為土質、地下水、氣候、坡角、坡面朝向、凍結層是邊坡失穩的主要影響因素。牛富俊等[11]指出，坡體表層土體變形和滑動可用有效應力分析法、總應力分析法、有效應力和固結理論等進行穩定性評價，認為地下冰水界面是熱融滑塌型滑坡的滑面，且基于靜力平衡法推導了安全系數表達式。

多年凍土區露天礦邊坡具有其特殊的地層分布，在水冰相變、凍脹融縮的周期性影響下易產生局部滑塌或掏蝕破壞。強度折減法、極限平衡法等多為整體穩定性分析評價方法，不適用于局部損傷破壞演化分析與評價。因此，本文擬沿用筆者前文提出的“含相變三場耦合簡化算法”和理想邊坡模型[12,13]，探究邊坡凍融局部損傷破壞演化過程，討論基于屈服接近度的局部穩定性分析評價方法的實用性，并以此研究各因素對邊坡穩定性的影響。

1 模型參數與方法

1.1 模型建立

木里煤田聚乎更礦區含F3、F7等斷層及其它軟弱結構面，但均埋藏較深，本文主要研究坡表局部的損傷破壞，故可忽略軟弱夾層等的影響。因此，沿用如圖2所示的理想邊坡模型[13]，采用FLAC3D軟件進行模擬分析。

圖2 計算模型Fig.2 Test model

其中，應力場、溫度場及滲流場分別采用M-C模型、各向同性傳導-對流模型、各向同性滲流模型。對模型底面進行全位移約束，側面進行法向位移約束。該地區坡體穩定，主要是坡表季凍層受溫度周期性的影響，而永凍層和不凍層相對比較穩定，基本不受大氣溫度影響，模型底面和側面與大地基巖的熱交換對坡表影響很小。為簡化起見，假設模型側面和底面與外界沒有流體和溫度交換，而坡表與外界有溫度交換，則坡表采用透水對流溫度邊界，其余采用不透水絕熱邊界。由邊坡巖土體自重產生初始應力場，并假設巖土體飽和。邊坡初始溫度T(h)按真實木里地溫且按深度(h)分段施加，坡表溫度T(d)按當地實際氣溫隨時間天(d)變化：

(1)

T(d)=-7.1-26.9×cos(2π(d+110)/12.0/30.0)

(2)

1.2 參數選取

凍巖試樣的相關試驗不少[2,12-19]，但該處邊坡涉及的是凍巖土體，其參數研究卻不多，因此只能假設、簡化或參見其它學者的研究成果。本文假設巖體參數隨溫度和凍融循環次數的變化規律近似同巖石試樣，則只須對初始值乘以對應的折減系數即可。假設巖體熱流場參數同巖石試樣，僅考慮溫度年凍融循環，而忽略溫度晝夜循環的影響。邊坡巖土體初始物理力學參數和熱流場參數分別如表1～2所示。各參數的取值，在相關文獻[12-13,16]中已進行了詳細討論。

表1 初始物理力學參數

表2 熱流學參數

凍融循環中巖石泊松比變化較小，可按初始值考慮[14]。假設只有變形模量和內黏聚力隨溫度和凍融循環次數[13]變化：

(3)

(4)

式中，T為溫度；n為凍融循環次數；E(T)、C(T)分別為溫度為T時對應的變形模量和粘聚力；KE、KC分別為凍融循環次數為n時對應的變形模量E和粘聚力C的凍融系數。

在整個凍融循環中，干燥巖樣表現為熱脹冷縮性，變形為彈性無殘余變形，但飽水巖樣經歷了冷縮-凍脹-回溫遲滯-融縮-熱脹五個階段，有明顯的殘余變形[15]。因此，對于飽水巖樣，巖基依然遵守熱脹冷縮規律，水分凍結體積增大使巖基孔隙擴大所產生的變形大于巖基自身冷縮變形，從而整體表現為凍脹；冰融化時體積變小，同時部分微裂隙閉合，從而整體表現為融縮；冰融化后巖基自身熱脹變形和不可恢復的蓬松化變形組成了其殘余變形。而對于凍巖土體工程而言，融化后水分散失，巖體強度也會降低，進而直觀表象是融沉或者滑塌破壞。飽水巖樣冰點上下分別表現“冷縮”和“冷脹”性質，且融縮后有殘余變形[15]。因此，巖石基質熱膨脹系數取3×10-6℃-1，水→冰相變時取-1.5×10-5℃-1，冰→水相變時可恢復熱膨脹系數為-1.2×10-5℃-1[12,14]。水、冰在冰點以上取2.1×10-4℃-1，在冰點以下取-1.5×10-3℃-1。由于溫度變化范圍相對較小，故可認為巖石基質的比熱和導熱系數不變，且對靜止空氣的對流傳熱系數可取5 W/(m2·℃)。

1.3 含相變THM耦合簡化算法

凍融循環下巖土體邊坡穩定是一個復雜的動態平衡問題，具有顯著的時效性，除了要考慮動態參數、邊界外，還須要考慮水冰相變、多場耦合的作用。本文沿用作者前文[12,16]中的“含相變THM耦合簡化算法”，對“熱-流-力”進行順次串聯得到耦合算法。其核心是水冰相變前后單元溫差對應的熱量累計不小于相變潛熱，只需采用FISH語言編寫簡單函數，將相變能量項用溫差能量項代換即可。即在水→冰相變時，系統溫度不變，釋放的潛熱抵消了系統在外界吸收熱量下引起的溫度降低；在冰→水相變時，系統溫度不變，儲存的潛熱抵消了系統在外界釋放熱量下引起的溫度升高。

1.4 屈服接近度分析方法

工程邊坡穩定性分析主要研究坡體破壞范圍、損傷破壞演化過程、變形破壞機理以及支護措施等，主要根據變形和應力判據來分析評價巖土體的安全狀態。而多年凍土區露天礦邊坡不同于一般邊坡，其具有特殊的地層分布，在季凍層底是相當厚度的永凍層，兩層之間則是富水的水冰臨界薄弱面。融化部分巖體強度降低，在自重作用下主要發生淺表層局部塊體或片狀脫落、滑塌、崩塌等破壞，可采用極限平衡法近似求解其表層局部穩定系數，也可間接利用邊坡局部經受凍融循環次數、局部變形量、局部損傷破壞體積等來評價，其中屈服接近度也是一種行之有效的方法。

周輝[19]提出屈服接近度(yield approach index，YAI)的概念，可廣義地表述為：描述一點的現時狀態與相對最安全狀態的參量的比。屈服接近度YYAI∈[0,1]，數值越小，表示越接近破壞，當YYAI=0表示應力點達到破壞。高麗燕[20]對YAI重新進行定義，空間應力狀態下一點的屈服接近度為：沿π平面坐標原點與應力點的連線方向，該點到屈服邊界線的距離與π平面坐標原點到屈服邊界線的距離之比，如圖3所示。

圖3 屈服接近度π平面示意圖Fig.3 YAI model on π plane

圖3中，AB為屈服邊界，Q點為O′P延長線與AB的交點，主應力空間的P點在π平面的坐標值可由這3個主應力計算得到。

對于M-C破壞準則，應力空間P點的屈服接近度可以表示為：

(5)

式中，σ1、σ2、σ3分別為P點第一、第二、第三主應力，上標′表示相應坐標軸名稱；σt為抗拉強度；Ι1為第一主應力不變量；J2為第二偏應力不變量；θσ為應力洛德角；c為內粘聚力；φ為內摩擦角。

2 邊坡穩定性分析

季節和晝夜交替產生的溫差引起巖體中的水分反復相變凍融，而反復凍融又嚴重削弱巖土體的強度，嚴重威脅邊坡的安全與穩定。多年凍土區露天礦邊坡穩定性影響因素不同于一般的邊坡，除了自身的形態尺寸外，主要還有溫度、水分、凍融次數、應力狀態等。因此，本節總結主要影響因素及其范圍，設計如表3所示的試驗方案，在基礎方案基礎上探究基于局部屈服接近度的穩定性分析評價方法的實用性及影響因素的敏感性。

表3 凍融循環作用下邊坡穩定性分析方案

其中，“邊界溫度升高量”和“凍結溫度降低量”是預設的兩個影響因素。全球變暖趨勢早有定論，且聯合國政府間氣候變化專門委員會(IPCC)估計，21 世紀全球平均氣溫將增加1.4～5.8 ℃，而青藏高原未來50a氣溫可能上升2.2～2.6 ℃[21]。

2.1 凍融破壞演化過程

對基礎方案的第5次凍融循環各月份的凍融深度、塑性區、拉應力區、剪應變增量區的演化過程進行分析，如圖4所示，由于受篇幅所限，僅列出1、4、7、10月的相關結果。

由圖4可知，3月底季凍層坡肩部位開始融化，7～9月達到最大融化深度約5～6 m，9月末氣溫下降進入負溫，季融層自上而下迅速回凍，至12月初季融層全部凍結，凍融深度與實際調研情況十分相似。坡表1～4 m范圍內主要存在拉剪塑性區歷史，有發生滑塌破壞的趨勢；坡表3～8 m范圍內主要存在拉塑性區歷史，可能由于水冰反復相變造成。回凍時由于凍脹力的作用，新拉塑性區會隨著凍結面的回遷而內移，而在坡肩部位新生成拉剪塑性區。當季凍層融化時，由于表層熱融滑塌的影響，會在坡肩、坡腳部位產生一定的拉應力。其中，8月拉應力區域最小，坡肩處拉應力不到0.05 MPa，而坡肩內一定深度卻達0.21 MPa；當季融層回凍時，拉應力區會逐漸擴大并向永凍層深入，但最大值部位由坡內向坡肩處外移，其中12月拉應力最大為0.50 MPa，則可近似認為坡肩處巖體的最大凍脹力為0.45 MPa。變形區主要在坡肩呈“三角形”狀分布，坡體淺層有沿季融層底冰面消融滑塌的趨勢，與頂蓋掏蝕現象類似，而坡肩內部是受力的關鍵，是坡體淺層安全的敏感區，如果在坡肩下部一定范圍內進行加固，可在很大程度上提高坡體安全性。

圖4 凍融演化過程Fig.4 Evolution process on freezing-thawing

2.2 局部屈服接近度

將YYAI平均分成10段(0.0～0.1,…,0.9～1.0)，在整個計算模型中，每段內各單元的總體積占模型總體積的百分比即分段體積百分比，將分段百分比逐段累加即基于屈服接近度分段體積累計百分比。對基礎方案的第5次凍融循環各月的屈服接近度分段的體積累計百分比進行分析，如圖5所示。

由圖5可知，各月份的屈服接近度分段體積累計百分比變化趨勢基本一致，而在YYAI=0.6的鄰域范圍內或曲線的轉折點附近各月份分布有所差異，其為易受水冰相變影響的區域，即凍融敏感區。且各月均在YYAI=0.8處增量最大，其中1月最小約為38%，而7月最大約為49%，其為坡體的內部穩定區。其中YYAI≤0.6的區域分布百分比約占整個模型體積的20%，基本屬于受外界環境影響的凍融區或季凍層，即有效屈服區。同時，高麗燕[20]指出，YYAI≥0.2表示結構安全，那么可認為YYAI≤0.6為凍巖土邊坡的損傷破壞的有效影響范圍，而其中0.2

圖5 屈服接近度分段的體積累計百分比Fig.5 Volume percentage on YAI distribution

雖然計算模型的邊界范圍取值對分段體積百分比有一定的影響，但對破壞損傷范圍的分布基本不影響，同時模型網格越小，范圍精度將越高。

YAI是空間應力狀態下定義的，是在網格單元尺度上判斷某一單元的損傷破壞程度，將同等級的YAI單元連接起來即YAI云圖或等值線圖。基礎方案的第5次凍融循環各月份(1、4、7、10月)的有效屈服接近度云圖，如圖6所示。

由圖6可知，有效屈服區YYAI≤0.6綜合了拉剪破壞性質及其區域分布，整體呈三角形分布，且破壞區YYAI≤0.2集中于坡肩部位，而損傷區基本沿凍融面成層分布，且隨著季融層的回凍在坡腳部位會產生屈服損傷區。

圖6 各月YAI演化過程Fig.6 Evolution process on YAI

YYAI表示的是某一單元的屈服或破壞程度，對于凍融邊坡工程的局部的安全程度而言，可采用邊坡表層易受凍融影響有效屈服區(YYAI≤0.6)內各單元的屈服接近度的體積加權平均值來表示，即局部屈服接近度：

(6)

其中，Yi和Vi分別表示各區段單元的屈服接近度和對應體積，且YJB值越低，表明損傷破壞程度越嚴重。對第5次凍融循環內各月份的局部屈服接近度進行分析，如圖7所示。

圖7 局部屈服接近度歷時曲線Fig.7 YJB changes with time

由圖7可知，隨著環境溫度的升高，季凍層逐漸融化，3月局部屈服接近度快速下降，4月至7月達到最低，而后隨著季融層回凍而升高。即每年4～7月份為危險月份，此時坡體淺層最可能發生局部崩塌或滑移破壞。因此，后續影響因素分析時，均以每年6月結果為基礎進行對比分析。由此可見，基于“局部屈服接近度”的方法來分析邊坡的局部安全穩定狀態可行，對解決此類邊坡問題有一定的參考意義。

2.3 影響因素的敏感性

為了方便后續影響因素分析，定義屈服相對比K，表示各對比方案的有效屈服范圍內破壞或損傷范圍分別與基礎方案的有效屈服范圍的比，K值越大，表明破壞或損傷程度相對越高。以第5次凍融循環6月的K值為例進行對比分析。

1) 影響因素整體分析

在基礎方案基礎上，在各影響因素變化范圍內，進行單因素分析，當各因素分別平均增加10%時，對應的K平均增量百分比，如表4所示。

表4 K平均增量百分比

注：表中加粗數值表示各區段的3個較大值。

由表4可知，當各因素分別平均增加10%時，凍結溫度對局部屈服區域影響最大，達到10.23%，而其中大部分主要影響的是損傷區，即永凍層上緣區域，那么在實際工程中應防止因人為因素而導致的凍結溫度降低。坡高和凍融循環次數對局部屈服區影響次之，且其依然主要影響的是損傷區，而坡角和熱膨脹系數直接導致坡體發生破壞。在參數變化范圍內，對局部有效屈服區的影響程度從大到小依次是：凍結溫度、坡高、凍融循環次數、熱膨脹系數、孔隙率，而其中對局部破壞區的影響程度從大到小依次是：坡角、熱膨脹系數、坡高、孔隙率、凍融循環次數。

2) 單因素影響分析

由上節可知，坡高、坡角、凍結溫度降低量、熱膨脹系數、孔隙率、凍融循環次數等對邊坡局部安全性影響較大，因此，本節對這些因素進行詳細分析，如圖8所示。

由圖8(a)、(b)可知，多年凍土深度范圍內，隨著坡高的增加，有效屈服區和損傷區逐漸增大，而已破壞和臨近破壞區先增加到40 m處達到極值而后逐漸降低。其中，坡高每增加10 m，已破壞、臨近破壞和損傷區比例分別增加0.009、0.017、0.093，說明坡高的增加容易導致坡體內部的損傷積累。隨著坡角的增加，有效屈服區變化不大，而損傷區逐漸減小，已破壞和臨近破壞區整體呈增大趨勢，而臨近破壞區先增加到50°處達到極值而后有降低趨勢。其中，坡角每增加10°，已破壞、臨近破壞和損傷區比例分別增加0.057、0.030、-0.093，說明坡角的增加容易導致坡體表層的破壞。

由圖8(c)、(d)可知，隨著凍結溫度的降低，有效屈服區和損傷區逐漸增大，而已破壞和臨近破壞區域變化不大。其中，凍結溫度每降低2℃，已破壞、臨近破壞和損傷區比例分別增加-0.001、0.010、0.332，說明凍結溫度的降低容易引起坡體內部的損傷積累。隨著熱膨脹系數的增加，各屈服區域均增加，其中熱膨脹系數每增加0.3 ℃-1，已破壞、臨近破壞和損傷區比例分別增加0.041、0.037、0.045，說明熱膨脹系數的增加將整體增大坡體的破壞與損傷區域。

由圖8(e)、(f)可知，隨著孔隙率的增大，有效屈服區、損傷區和已破壞區先增加到0.2處達到極值而后逐漸降低。其中，在極值前孔隙率每增加0.05，已破壞、臨近破壞和損傷區比例分別增加0.050、0.026、0.043，說明孔隙率的增加將整體增大坡體的破壞與損傷區域。隨著凍融循環次數的增加，各屈服區域均呈增加趨勢但增率變緩，且早期主要影響的是已破壞區域而后期主要影響損傷區域。其中，凍融次數每增加5次，已破壞、臨近破壞和損傷區比例分別增加0.011、0.014、0.136，說明凍融次數的增加主要引起損傷區的累積。

3 結 論

本文沿用 “含相變三場耦合”簡化算法和理想邊坡模型，探究了邊坡凍融局部損傷破壞演化過程，討論了基于屈服接近度局部穩定性分析方法的實用性，并以此探究了各因素對邊坡穩定性的影響。

1) 模擬的凍融深度與實際調研情況十分相似，新拉塑性區會隨著季融層回凍而內移，而拉應力最大值區域反而外移至坡肩。變形區主要在坡肩呈“三角形”狀分布，有產生消融滑塌的趨勢，與頂蓋掏蝕現象類似。坡體淺層有沿季融層底冰面滑塌的趨勢，而坡肩內部是坡體整體安全的受力關鍵。

2) 凍巖土體邊坡的局部有效屈服區為YYAI≤0.6，而其中0.2

3) 各影響因素對局部有效屈服區的影響程度從大到小依次是：凍結溫度、坡高、凍融循環次數、熱膨脹系數、孔隙率，且坡高、坡角、孔隙率分別存在40 m、50°和 0.2的極值。坡高的增加、凍結溫度的降低、凍融次數的增加主要導致坡體內部的損傷積累，坡角的增大更易導致坡體表層的破壞，而熱膨脹系數的增加、孔隙率的增加將整體增大坡體的破壞與損傷區域。