水-熱耦合下露天礦凍結期靠幫開采邊坡穩定性研究

谷浩源 ，韓 流 ，張 健 ，劉世寶 ，高志強 ，吳鋒鋒

（中國礦業大學 礦業工程學院，江蘇 徐州 221116）

煤炭作為我國的主要能源，其開采量和消費量均居能源結構的首位[1]。我國目前擁有露天礦山20 000 余座，超過80%的露天煤礦位于北方季節性凍土地區，凍結期大多在半年以上[2-3]。以內蒙古自治區的露天煤礦為例，因其受季節性因素影響較大，其礦區生產在冬季和夏季對應剝離作業期和剝離停產期，且露天礦邊坡處在反復凍融條件下，邊坡穩定性也受一定的影響。為了提高采出率、降低剝采比，實現煤炭資源的回收和經濟效益的提高，靠幫開采技術得到了廣泛的應用[4]。冬季靠幫開采回收端幫壓煤，可通過內排壓幫提高幫坡角度，幫坡角的變化直接影響煤炭采出率，而邊坡穩定性是保證安全開采的前提，因此這兩者是靠幫開采的關鍵參數。

針對邊坡穩定性和幫坡角的研究，國內外許多學者進行了研究，取得了豐碩成果。吳榕真等[5]利用FLAC3D軟件，結合不同靠幫角度，分析了靠幫開采過程中邊坡的穩定性變化，確定了靠幫開采最佳的端幫邊坡角；佘長超等[6]對排土場邊坡軟巖物料進行了凍融循環下剪切破壞試驗，得出了凍融循環次數增加將導致邊坡穩定性下降的結論；馬忠輝等[7]依據冬夏兩季邊坡巖層的物理力學參數，運用穩定性分析軟件，計算兩季端幫穩定性系數強化幅度，確定出不同季節下靠幫開采方案；陸翔[8]結合寒區季節性氣候變化特征，分析了凍融作用對端幫邊坡穩定系數及破壞模式的影響，基于凍結泥巖強度強化特性，提出了寒區露天礦季節性靠幫開采方法。雖然上述學者對露天煤礦端幫靠幫開采進行了深入研究，但針對凍結期凍結層厚度變化對靠幫開采邊坡穩定性以及幫坡角的影響，尚未形成統一的認知。

基于此，以寶日希勒露天礦北部端幫為工程背景，基于多孔介質理論和水-熱耦合理論，考慮凍結期不同厚度的凍結層對邊坡穩定性影響的變化規律，建立軟巖邊坡穩定性數值計算模型，揭示凍結期靠幫開采邊坡穩定性變化規律以及合理的開采時間。

1 水熱耦合模型構建

隨著露天礦進入凍結期，溫度的改變將會引起巖土體內部冰水相變以及水分遷移，而冰水相變所產生的潛熱會反作用于溫度場，且相變引起的變形也將使得土體力學參數發生變化[9]。針對水-熱耦合兩相物理場問題，可通過構建二者獨立方程，尋求耦合聯系方程，為凍結期露天礦邊坡水熱特性分析提供理論依據。

1.1 溫度場控制方程

基于傅里葉熱傳導定律和能量守恒原理，考慮二維水熱耦合問題，并考慮相變熱，針對多孔介質凍巖土導熱控制微分方程如下[10]：

式中： ρ為巖、土密度，kg/m3；C(θ)為巖、土體積比熱容，J/（kg·℃）[2,11]；T為巖、土瞬態溫度，℃；t為時間，s； λ(θ)為巖、土導熱系數，W/（m·℃）；?為微分算子，對于二維問題為[?/?x,?/?y]；L為相變潛熱，取值為334.5 kJ/kg；ρi為冰的密度，kg/m3； θ為凍巖、土體積含水率；θi為凍巖、土中冰的體積含量。

1.2 滲流場控制方程

邊坡巖土凍融過程中，始終存在自由水流動，其運動過程符合廣義達西定律。考慮邊坡凍融水分場運動為非飽和多孔介質流動，選用Richards水分遷移方程[12]。考慮冰水相變對滲流場的影響，得凍土中未凍水遷移運動方程：

式中： θu為凍巖、土中未凍水的體積含量；ρw為凍土中水的密度，kg/m3；D(θu)為季凍邊坡巖、土體的水力擴散系數；k(θu)為季凍邊坡巖、土體在重力加速度方向的土體滲透系數。

為充分描述季凍巖、土體邊坡滲流真實情況，TAYLOR[13]在確定水分擴散系數時，引入阻抗因子I表示季凍巖、土體孔隙中冰對未凍水遷移的阻礙，并給出了水分擴散系數的公式：

為簡化計算引入相對飽和度S：

式中：c(θu) 為 比水容量，1/m；kS為飽和巖土滲透系數，cm/s；a0、m、l為與巖、土本構相關的參數，取值參考Van Genuehten（VG）滯水模型[14]；θr為 殘余含水率； θs為飽和含水率。

固液比BI相關概念的提出[11]，使得季凍邊坡巖、土體由冰水相變引起的溫度、體積含水量和體積含冰量三者變化產生關聯。其計算式如下：

式中：Tf為初始土體凍結溫度，℃；B為土質含鹽量變化經驗系數[15]。

1.3 模型及驗證

基于COMSOL Multi-physics 有限元數值模擬軟件，具有的多物理場耦合能力以及非線性微分方程組求解能力，通過對其進行二次開發，實現邊坡巖、土體水熱力場耦合模擬，軟件中內置系數型偏微分方程（PDE）形式如下式：

式中：ea為方程質量系數；u為方程自變量；da為 阻尼或質量系數；c為擴散系數； α為守恒通量對流系數； γ為守恒通量源； β為對流系數；a為吸收系數；f為源項。

將式（7）、式（8）代入式（1）中，整理后得溫度場控制方程為：

為驗證模型開發的合理性，參考魏道凱[16]對一維土柱凍脹試驗的數據進行測試。土柱寬0.14 m，高0.4 m，采用與原文相同的四邊形單元進行網格劃分，水分場溫度場邊界條件與原文一致。模型計算所需參數為[17]：①材料：黏性土；②Tf：-1.33 ℃；③B：0.56；④a0：2.65；⑤m：0.25；⑥l：0.5；⑦S：0.76；⑧θs：0.432；⑨θu：0.33；⑩ks：10-8cm/s。

不同時刻溫度和未凍水含量的試驗值與模擬值結果對比如圖1。

圖1 模型試驗值與模擬值對比Fig.1 Comparison of model test values and simulation values

從結果來看，模擬值得出的溫度和未凍水含量變化規律與試驗值一致，都呈現隨時間增加而下降的趨勢，而且無論是溫度還是未凍水含量變化，在同一時刻其試驗值均大于模擬值，這是由于數值模擬相較于實際試驗而言，過于理想化，實際試驗無法做到保溫隔熱，避免能量損失。綜上所述，所采用的水熱耦合計算模型與實際試驗數據整體上貼切，證實了所開發的水熱耦合模型描述凍土水熱特性的合理性。

2 凍結期邊坡穩定性

北方露天礦多以土質邊坡為主，每年10 月份開始進入凍結期，巖土孔隙水受溫度影響冷結為冰。隨著低溫時間延長，凍結層向深部擴展形成季節性凍土。此過程主要分為初冬凍融交替、凍結發展以及融凍3 個階段[18]。研究凍結期（凍結發展階段）凍結層厚度對邊坡穩定性的影響，可為凍結期露天礦靠幫開采邊坡穩定性分析提供理論依據。

2.1 凍結期巖土體強度參數變化

土和巖石都屬于多孔介質。以土為例，土是由固體顆粒、水和氣體組成。非凍結期和凍結期土體結構如圖2。

圖2 非凍結期和凍結期土體結構Fig.2 Soil structure in non-freezing period and freezing period

天然狀態下，土是三相體系，而在形成凍土后就變為“四相”體系。凍土形成過程中，溫度的改變將會引起冰水相變以及水分遷移，冰水相變所產生的潛熱會反作用于溫度場，相變引起的變形也將影響土體力學參數，而且凍結后土的強度增加往往是大于其天然狀態下強度和冰強度之和[19]。對于土質邊坡穩定性而言，抗剪強度是主要的影響因素，其中黏聚力和內摩擦角是關鍵參數[20]。寶日希勒露天礦區凍結期和非凍結期地層物理力學參數見表1。

表1 非凍結期和凍結期邊坡地層物理力學參數Table 1 Physical and mechanical parameters of slope strata in non-freezing period and freezing period

2.2 凍結期溫度對凍結層厚度的影響

根據寶日希勒露天礦采掘平面圖選擇某一北部端幫剖面，通過在Auto-CAD 中畫出寶礦北幫凍結期模型，模型長600 m，高253 m，進入凍結期后，邊坡臺階會因為外部溫度降低的影響形成凍土層，凍結期北幫模型示意圖如圖3。

圖3 凍結期北幫模型示意圖Fig.3 North slope model of freezing period

研究凍巖土試樣的試驗不少，但涉及邊坡凍巖土體的研究不多，研究參數相對而言較為稀缺，因此只能假設、簡化或參見其它學者的研究成果。假設邊坡巖體熱力學和水力學參數同巖石試樣，部分巖體熱、水力學參數可用物理力學參數近似的巖樣代替。各參數的取值參考相關文獻中進行了詳細討論[12,21-24]，北幫邊坡水熱學參數見表2，北幫邊坡水力學參數見表3。

表2 北幫邊坡水熱學參數Table 2 Hydrothermal parameters of north slope

表3 北幫邊坡水力學參數Table 3 Hydraulic parameters of north slope

模型整體采用三角形單元進行網格劃分，對局部加密處理，保證計算結果的精確性。根據地質資料，設置初始相對飽和度為0.6，模型水分場邊界條件設為0 通量，兩側溫度邊界條件為絕熱邊界，下邊界為溫度與當地年平均氣溫保持一致的恒溫帶，模型上邊界采用Dirichlet 溫度邊界條件，參考寶日希勒露天礦當地年平均氣溫，選擇10 月中旬作為凍結期起點，得出上邊界溫度T（t）表達式：

式中：t為凍結期時間，d。

計算結果表明，坡面各處凍結厚度大致相同，以最低凍結溫度-0.54 ℃所在位置作為凍結層厚度判斷標準，最大凍結厚度不超過5 m。為此，選擇高程+252 m 坡頂面中點為特征點，沿坡面內法線方向做向下延伸5 m 的截線，提取不同時間截線的溫度作為凍結層厚度的判斷依據。不同時期的邊坡凍結層厚度如圖4。

圖4 不同時期的邊坡凍結層厚度Fig.4 Thickness of frozen layer of slope at different periods

由圖4 結果可知：初始凍結期為0 d 時，即每年10 月中旬，此時邊坡整體溫度為5 ℃，未產生凍結；隨著時間增加，在凍結期60 d，即12 月中旬時，邊坡凍結層厚度從0 m 增加到1.82 m；在前145 d 內，凍結層厚度隨時間不斷增加屬于凍結發展階段，最大可達3.45 m，這與寶日希勒露天礦地質勘探報告中最大凍結厚度3.5 m 接近，印證了模型開發的合理性；從145 d 至160 d，凍結層厚度維持在3.4 m 左右，出現穩定凍結；隨著氣溫的逐漸回升，180 d 時邊坡整體凍結厚度出現不連貫現象，在厚度為2.38 m 到3.13 m 處溫度達到了凍結溫度，凍結層包含在邊坡內部。這是因為此時邊坡處于融凍階段，初期凍土層停止向下發展，凍土層厚度達到穩定，后期凍土層由底部和表層向凍土層消融解凍，在坡體內部形成不連貫凍結層，因此可以確定160 d 以后邊坡處于融凍階段，直至190 d 后，即第2 年4 月中下旬邊坡凍層完全融化。

2.3 不同凍結層厚度下邊坡穩定性

目前，關于邊坡穩定性的計算方法，主要有極限平衡法、概率分析法、數值計算法[25]。其中，極限平衡法是工程領域中最基本的方法。GEOslope 是工程領域用于分析邊坡穩定性十分常用的軟件[26]，其主要是依據極限平衡法的分析原理，對邊坡穩定性系數進行求解。基于穩定性系數Fs的極限平衡法有多種，常用于邊坡穩定性分析的 方 法 有Morgenstern-Price 法[27]、Spencer 法[28]、JanbuGeneralized 法[29]、Bishop 法[30]。因礦山實際生產情況復雜，選擇上述4 種方法進行計算，根據不同方法結果進行對比，選取較為合適且相對保守的穩定系數作為邊坡穩定系數。

利用水熱耦合模型得出的不同時期凍結層厚度，結合GEO-slope 軟件進行凍結期邊坡穩定系分析。利用不同方法對北幫邊坡進行不同時期不同厚度凍結層下邊坡穩定性分析，不同凍結層厚度下北幫穩定性系數見表4。

表4 不同凍結層厚度下北幫穩定性系數Table 4 Stability coefficients of north slope with different thickness of frozen layer

由表4 可知：隨著凍結時間的增加，凍結層的發育，無論是哪種計算方法，其邊坡穩定性都是在增加的，說明在凍結期巖體強度的提升，凍結層厚度的變化將對邊坡穩定性有一定影響。依照選取較為合適且相對保守的穩定系數為選取原則，可以得出凍結期靠幫開采最大穩定系數為1.336，相較于非凍結期1.245 提升了7.3%。從表4 中可以看出邊坡穩定性系數在2 月中旬到3 月中旬增幅為0.07%，可認為在此期間的穩定性系數最大且變化不大，靠幫開采的時機最佳。

3 凍結期靠幫開采幫坡角度

韓流等[31]、尚濤等[32]針對凍結期靠幫開采提出了條分式靠幫開采方法，在保證邊坡穩定的前提下，提高了煤炭回收量的同時降低了剝采比，結合寶日希勒露天礦北幫實際情況凍結期選擇不剝離覆巖的靠幫開采[20,33]。

寶日希勒礦區每年的10 月初進入凍結期，于第2 年4 月中旬解凍，最大凍土厚度可達3.45 m。選擇在凍結2 月中旬至3 月中旬期間進行靠幫開采，此時的邊坡穩定系數最高，靠幫開采角度相對而言較大，端幫煤炭回收量較多。通常選取2倍工作幫臺階高度作為內排壓幫的臺階高度，采用條分式靠幫開采，選擇30 m 作為壓端幫底腳臺階的寬度，并用與端幫底腳高度相同的內排臺階及時跟進壓幫，臺階高度為25 m。構建的端幫穩定性分析模型如圖5。

圖5 凍結期北幫穩定性分析模型Fig.5 Analysis model of north slope stability during freezing period

寶日希勒露天礦原北幫幫坡角為26°，當凍結層厚度為3.4 m 時穩定系數為1.336。分別對幫坡角為27°、28°、29°以及30°的端幫進行穩定性計算并得出穩定系數，得出的北幫在不同幫坡角時邊坡穩定性系數見表5。由表5 可知：隨著幫坡角不斷增大，其對應的穩定系數會隨之減小。

表5 不同幫坡角的穩定系數Table 5 Stability coefficients of different slope angles

根據邊坡穩定性的判定依據可知，當滑體上的抗滑力（矩）與下滑力（矩）比值大于1，即抗滑力大于下滑力，此時滑體處于安全狀態，邊坡穩定；反之，邊坡失穩。在滿足規范端幫邊坡穩定系數1.10 前提下，確定靠幫開采的最佳幫坡角為29°。

4 工程實踐

根據上述研究可知，凍結期內由于邊坡表面形成了強度較大的凍結層，邊坡整體穩定性提升，穩定系系數得到增大。在此基礎上，對寶日希勒露天礦凍結期采用內排壓幫寬度30 m，內排臺階高度25 m 的不剝離靠幫開采方式，進行凍結期靠幫開采。靠幫開采合理幫坡角為29°，通過CAD計算得可回收的煤層橫截面面積為180 m2，根據地質資料描述，采用350 m 的可采長度計算，可回收煤炭8.4 萬t，以 150 元/t 的市場價格計算，可創造1 260 萬元的經濟價值，由此可見基于凍結期軟巖強度提升的靠幫開采產生的經濟效益是十分顯著的。

融凍階段坡體含水量如圖6。

圖6 融凍階段坡體含水量Fig.6 Water content of slope in thawing and freezing stage

冬春交替時節，氣溫逐漸上升，當最低氣溫達到0 ℃及以上，露天礦邊坡巖土進入解凍期，如圖6（a）；隨著時間增加，氣溫回暖，凍結層開始融化使得邊坡巖體含水量增大，坡體水壓力增大，滑坡危險系數增加，如圖6（b）。

解凍層巖石黏聚力為20.5 kPa、內摩擦角為32.3°，密度為2.20 t/m3。將解凍層巖石力學參數代入到Geo-slope 中，得出解凍期幫坡角為29°時端幫穩定系數為1.084，無法達到規范端幫邊坡穩定系數1.10 的要求。為滿足安全要求，可通過增加內排臺階寬度的方法，使其內排壓幫體積增大，以保證解凍后端幫穩定系數的增加。解凍期端幫保障措施如圖7，圖中△V為增大內排臺階寬度所需的廢棄物量。

圖7 解凍期端幫保障措施Fig.7 Guarantee measures for north slope in thawing period

臺階寬度以5 m 為1 個步距遞增，經多次穩定性計算后，確定內盤寬度在原有基礎上增加10 m，即內排寬度為40 m 時，此時的穩定系數為1.122，滿足規范端幫邊坡穩定系數1.1 的要求。

5 結 語

1）針對凍結期凍結層厚度受溫度、水分變化影響，利用水熱耦合控制方程構建了1 個分析凍結期邊坡凍結層厚度變化的模型。

2）邊坡穩定性系數與凍結層厚度呈現正相關，凍結期內穩定性系數隨著凍結厚度增加而增大，降低而減小。靠幫開采最佳時機為2 月中旬至3月中旬，凍結層厚度達到3 m 以上，最大為3.45 m，穩定性系數1.336。

3）得出了最佳時機進行靠幫開采，合理幫坡角為29°，與原設計26°相比，多回收煤炭8.4 萬t，創造1 260 萬元經濟價值，考慮解凍期邊坡穩定性降低，通過增大內排臺階寬度40 m 的方式，保障了邊坡過渡期安全。