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

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

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

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

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

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

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

多年凍土區(qū)露天礦邊坡具有其特殊的地層分布，在水冰相變、凍脹融縮的周期性影響下易產(chǎn)生局部滑塌或掏蝕破壞。強度折減法、極限平衡法等多為整體穩(wěn)定性分析評價方法，不適用于局部損傷破壞演化分析與評價。因此，本文擬沿用筆者前文提出的“含相變?nèi)龍鲴詈虾喕惴ā焙屠硐脒吰履Ｐ蚚12,13]，探究邊坡凍融局部損傷破壞演化過程，討論基于屈服接近度的局部穩(wěn)定性分析評價方法的實用性，并以此研究各因素對邊坡穩(wěn)定性的影響。

1 模型參數(shù)與方法

1.1 模型建立

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

圖2 計算模型Fig.2 Test model

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

(1)

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

(2)

1.2 參數(shù)選取

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

表1 初始物理力學參數(shù)

表2 熱流學參數(shù)

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

(3)

(4)

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

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

1.3 含相變THM耦合簡化算法

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

1.4 屈服接近度分析方法

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

周輝[19]提出屈服接近度(yield approach index，YAI)的概念，可廣義地表述為：描述一點的現(xiàn)時狀態(tài)與相對最安全狀態(tài)的參量的比。屈服接近度YYAI∈[0,1]，數(shù)值越小，表示越接近破壞，當YYAI=0表示應力點達到破壞。高麗燕[20]對YAI重新進行定義，空間應力狀態(tài)下一點的屈服接近度為：沿π平面坐標原點與應力點的連線方向，該點到屈服邊界線的距離與π平面坐標原點到屈服邊界線的距離之比，如圖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為內(nèi)粘聚力；φ為內(nèi)摩擦角。

2 邊坡穩(wěn)定性分析

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

表3 凍融循環(huán)作用下邊坡穩(wěn)定性分析方案

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

2.1 凍融破壞演化過程

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

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

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

2.2 局部屈服接近度

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

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

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

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

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

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

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

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

(6)

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

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

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

2.3 影響因素的敏感性

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

1) 影響因素整體分析

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

表4 K平均增量百分比

注：表中加粗數(shù)值表示各區(qū)段的3個較大值。

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

2) 單因素影響分析

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

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

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

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

3 結(jié) 論

本文沿用 “含相變?nèi)龍鲴詈稀焙喕惴ê屠硐脒吰履Ｐ?，探究了邊坡凍融局部損傷破壞演化過程，討論了基于屈服接近度局部穩(wěn)定性分析方法的實用性，并以此探究了各因素對邊坡穩(wěn)定性的影響。

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

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

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