高原多年含水凍巖爆破炮孔溫度-水分場耦合分析

費鴻祿，甄 帥，李文焱，胡 剛，聶 寒

(遼寧工程技術大學 爆破技術研究院,阜新 123000)

高原高寒地區的露天礦山進行爆破作業時，爆破區域炮孔成孔后，滲流水通過巖石裂隙流入孔內，由于受當地海拔高度、氣候、晝夜溫差、水文地質的多種因素影響，經過一定時間導致炮孔內滲流水完全凍結，從而無法進行裝藥工序。研究表明，影響滲流水凍結的主要因素分別為場區溫度場作用以及滲流水速度，相比而言，溫度場所引起的效應更大[1,2]。

國內外許多學者對凍土溫度場已經有了較為深入的研究，Noorishad等基于Biot固結理論[3]，將變分原理和伽遼金公式與有限元方法相結合，首次提出了飽和巖體固-液-熱耦合方程。此后，胡向東等通過單排管、雙排管、三排管對凍結壁溫度影響的一系列研究[4-9]，提出“卸載狀態下凍結壁-周圍土體共同作用”概念并建立相關凍結壁力學模型，推導出凍結穩態溫度場解析解矩陣表達式。王效賓等以南京地鐵線工程為例[10]，采用ADINA有限元軟件分析地鐵線溫度場，得出了導熱系數、比熱容等因素對凍土溫度場的影響。張松等為研究突發滲流作用下溫度場溫度演變過程基于相似模型開展凍結試驗并得到相關規律[11,12]。邵玉龍等建立了描述三維裂隙網絡巖體滲流傳熱耦合模型[13]，并與二維裂隙傳熱模型做對比，發現大量流體流經基巖時會通過裂隙同周圍巖體進行熱量交換，造成裂隙周圍巖體溫度首先發生改變，然后向周圍擴散的規律。黃峰等通過有限元軟件ANSYS和TAITHERM建立了武漢某立交橋的三維溫度場計算模型[14-16]，數值模擬結果與實時監測數據耦合度頗高，為類似工程溫度場研究提供參考依據。邰博文等通過COMSOL有限元模擬及現場試驗揭示了寒區鐵路路基的水熱耦合機理[17-20]，并分析其與時空演變規律。韓小妹等在高寒地區通過溫控監測數據和研究分析成果[21]，總結出混凝土壩溫度場和應力場的變化規律。周曉敏等在礦井工程中采用數值模擬和現場試驗對滲流作用下凍結壁的溫度場發展規律進行研究[22,23]。Yuchao Zeng等研究了地熱田溫度場、地下水初始靜壓力場和水密度分布特征[24]，同時分析影響溫度場的主要因素，并考慮了裂縫系統的地質特征，建立了裂縫系統的概念模型。

綜上所述，目前國內外學者在分析溫度場-水分場時，大多是基于ANSYS、ADINA等有限元軟件或現場數據及理論公式進行單一變量分析并得出相應結論，并且主要針對隧道凍結壁、深井工程、大體積混凝土等大深度工程進行溫度場研究，對高寒常年凍巖地區露天礦山爆破作業過程中解決炮孔出現涌水、冰凍現象對溫度-水分場多場耦合做深入系統研究的鮮見。因此，本文以拉薩市墨竹工卡縣甲瑪礦為背景，以光纖測溫儀的試驗結果與COMSOL-Multiphysics有限元模擬相結合，對比分析高原地區炮孔成孔后溫度場發展規律。通過研究溫度場的時空關系，從根本上解決常年凍巖條件下炮孔冰凍裝藥的技術難題，為高原露天礦高效開采提供有效的技術方法。

1 工程背景

甲瑪露天礦位于西藏拉薩市墨竹工卡縣，年度平均氣溫2～17℃，高寒干燥，晝夜溫差大，在高寒高海拔條件下，最低溫度可達-10℃。在進行爆破作業時，導致爆破場地內炮孔常年出現地表向孔內滲水和鉆孔時孔壁涌(滲)水、冰凍現象。主要表現為：在炮孔成孔后3～12 h孔內滲流水經過礦區溫度場的重分布在短時間內可完全凍結，導致大量完成的炮孔無法實施裝藥工序，當前應用最廣泛的解決措施是選擇傳統的鉆機透孔法進行破冰處理，但此方法在實際應用中存在多個不足，首先此方法是通過壓縮透孔時間以滿足緊后工序的裝藥步驟，這樣就有可能出現交叉作業的風險增加，造成安全隱患。其次二次透孔的成本也是不可忽略的，增加透孔工序會嚴重制約爆破作業的效率。

因此研究炮孔成孔后溫度場變化規律，并根據研究成果開發相對應的技術、裝備、方法，對爆破工序的正常進行具有重要意義。為研究高原多年含水凍巖溫度場的發展規律，以拉薩市墨竹工卡縣甲瑪露天礦5120平臺為工程試驗背景，如圖1所示。

圖 1 甲瑪露天礦Fig. 1 Jiama open-pit mine

2 耦合理論分析

2.1 溫度場控制方程

溫度場是指某研究物體在某時間內所處空間溫度的集合，它反映了溫度的時空間關系，通常用空間直角坐標系以及時間變量表示,其方程如下所示

T=f(x,y,z,t)

(1)

式中：T表示溫度；x、y、z為空間直角坐標系坐標；t為時間。

礦區炮孔內溫度場為三維熱傳導問題，但由于炮孔成孔后為空心圓柱，因此可取對稱結構溫度變化，簡化為二維熱傳導問題。如(2)所示

T=f(x,z,t)

(2)

假設巖體為各向同性且均勻的，通過數學物理方法中熱傳導方程拉普拉斯公式，基于彈性力學理論推導出其二維熱傳導微分方程

(3)

(4)

式中：α為常數；Q表示單位熱量，J；c表示物體比熱容，J/(kg·K)；ρ表示物體密度，kg/m3；2為拉普拉斯算子。

考慮相變潛熱的傳導方程

(5)

式中：C為熱容量，J/(kg·K)；θ為體積含水量，kg/m3；λ為導熱系數，W/(m·K)；對于簡化二維問題為為相變潛熱；ρI為冰的密度，kg/m3；θI為冰的體積含量，kg/m3。

在進行有限元軟件COMSOL-Multiphysics數值建模時，炮孔內滲流水由液相轉變為固相需要放熱，其釋放的熱量會通過巖石孔壁及孔口空氣流動發生熱傳導。

2.2 水分場控制方程

在炮孔成孔后礦區原地質構造破壞會產生裂隙，裂隙中存在的水通過裂隙滲流至炮孔，其滲流作用符合達西定律。根據Richard方程并考慮冰與水的相變可得非飽和巖石內水分場遷移方程

(6)

式中：ρW表示水的密度，kg/m3；θu為液態水體積量，kg/m3；k為滲透系數。冰凍巖土中水的擴散率如下

(7)

式中：k(θu)為土體滲透率，m/d；c(θu)為比水容量，kg；I為阻抗因子[25]，表示孔隙中冰對水產生滲流的阻礙作用。

2.3 溫度場與水分場的耦合

溫度場控制方程中，溫度T為自變量。水分場控制方程中，液態水體積含量θu與冰體積含量θI為自變量。在有限元軟件建模求解時，兩個控制方程不足以將方程中三個未知變量表示，需要引入耦合項將水熱方程聯系求解，因此選用“固液比”的概念作為耦合項[26]。固液比即為巖土中冰的體積含量與自由水體積含量的比值，其公式可表示為

(8)

式中：Tf表示巖土體凍結溫度，K；B為常數。

因此冰的體積含量為

θI=BI·θu

(9)

由巖土體的相對飽和度公式，結合VG滯水模型及Gardner滲透模型[27]，可推導出水的體積含量

(10)

θu=(θs-θr)·S+θr

(11)

式中：s為相對飽和度；θs為飽和含水率；θr為殘余含水率。因此溫度場控制方程中熱源項可化簡為

(12)

引用“固液比”的概念后，溫度場控制方程為

(13)

水分場控制方程為

[D(S)S+k(S)]

(14)

此時可基于COMSOL-Multiphysics有限元軟件建立溫度-水分耦合模型。

3 數值模型建立

3.1 選擇物理場

為分析礦區炮孔涌水凍結過程及炮孔內溫度場發展規律，考慮到溫度場-水分場的相互耦合作用，在COMSOL-Multiphysics的模型庫中選擇多孔介質傳熱模塊，并添加流體域及固體域進行有限元分析。為以下敘述方便，稱炮孔成孔后至裝藥前之間的時間段稱之為空孔期，結合本工程背景的實際，空孔期孔內溫度發展規律對于涌水凍結具有決定性影響作用，現場實際空孔期在24 h左右，因此對溫度場進行裝藥時間最大間隔24 h內的瞬態研究。

3.2 建立幾何模型

使用COMSOL-Multiphysics組件中幾何模型建立板塊。在空孔期中孔內溫度主要受到自由水在炮孔壁裂隙的滲流作用影響，同時，炮孔內的溫度場通過空氣對流以及孔內與地表面的熱交換效應，也會對炮孔內的溫度變化造成一定影響。簡化炮孔成孔后孔內y方向溫度變化，建立相應二維幾何模型，模型總分為2部分，中間部分為直徑0.12 m、深度18 m炮孔，炮孔周邊為待爆礦區，左右區域各設5 m影響范圍，如圖2所示。

圖 2 礦區幾何模型Fig. 2 The geometric model of the mining area

3.3 模型邊界條件的確定

(1)力學邊界條件的確定

模型兩邊為輥支撐，底面為固定約束，頂面為自由約束。

(2)溫度邊界條件的確定

炮孔成孔后溫度主要受炮孔周邊巖體的熱傳導作用影響，使用SG-DTS-84U光纖測溫儀進行礦區穩定狀態下孔內溫度測量作為炮孔成孔后的初始溫度。取不同地表溫度分析初始地表溫度對炮孔溫度發展的影響，并以外部自然對流的方式向礦區底部傳熱，礦區底面平均溫度為炮孔底部溫度。

(3)水分場邊界條件的確定

炮孔在成孔后裂隙水會隨著孔壁流入孔內，由于滲流速度在不同平臺下各不相同，因此分析不同滲流速度下裂隙水對炮孔內溫度場的影響，并給予流體添加相變材料，由液態轉變為固態[28]，相變溫度273.15[K],相變轉變間隔2[K]，相變潛熱為333[kJ/kg]。見圖3。

圖 3 相變潛熱Fig. 3 Latent heat of phase change

3.4 模型材料定義

模型材料分為流體、固體、冰3種。具體屬性如表1所示。

表 1 模型材料

3.5 劃分網格

使用COMSOL-Multiphysics內置Mesh對模型網格劃分，全局劃分定義超細化，礦區邊界部分定義為極細化，炮孔內部定制單元大小參數，最大面積單元0.06 m2，最小面積單元大小0.0004 m2?？傆嬎銌卧?6392個，計算網格面積202.4 m2。如圖4所示。

圖 4 幾何模型網格劃分Fig. 4 Meshing of geometric model

4 模擬結果與分析

拉薩市墨竹工卡縣冬季平均氣溫為-3～13℃，夏季平均氣溫6～20℃，由于甲瑪礦區處在高寒高海拔地區，因此溫度相對于墨竹工卡縣更低?？紤]礦區地表區域溫度受外部環境影響大，炮孔成孔后滲流水的流速對溫度的擴散起決定性作用，因此本模型分別對外界溫度為0℃(273.15 K)、5℃(278.15 K)、10℃(283.15 K)滲流水在炮孔內的涌出速度1 m/h、2 m/h、3 m/h、5 m/h、10 m/h的情況進行模擬分析，當滲流水溢滿炮孔時水位達到相對穩定狀態不再上升。使用COMSOL-Multiphysics有限元分別分析成孔后礦區溫度場變化情況，并導出各個時間段的溫度數據。

4.1 炮孔成孔后24 h內溫度分析

結合現場工程情況，對地表溫度為5℃(278.15 K)炮孔內滲流速度為3 m/h時成孔后24 h內溫度發展規律進行模擬分析。t=0時礦區溫度分布如圖5所示，滲流水速度3 m/h時炮孔成孔后1 h、2 h、4 h、8 h、12 h、16 h、20 h、24 h孔內溫度場隨時間分布圖如圖6所示。通過有限元軟件模擬結果，對比分析炮孔成孔后24 h內溫度隨深度變化規律，如圖7所示。

圖 5 成孔后礦區初始溫度分布圖Fig. 5 Initial temperature distribution in the mining area after blast hole formation

圖 6 滲流水速度3 m/h溫度隨時間變化圖Fig. 6 The temperature changes with time with the seepage water velocity of 3 m/h

由圖5可見，炮孔成孔后t=0時礦區溫度分布均勻，且溫度隨深度遞減。

由圖6模擬結果顯示：在炮孔成孔后，孔內溫度與礦區內部溫度隨時間變化明顯，炮孔內0～8 h溫度變化快，8 h后逐步趨于穩定，其相同深度的溫度差異不大。炮孔深度0～5 m區域溫度受對流熱通量影響大，在深度3～6 m時呈現出溫度分層。深度5～18 m區域溫度變化由炮孔周邊巖體溫度決定，溫度隨深度整體呈下降趨勢。

由圖7可見：成孔1 h內與炮孔深度5 m左右位置存在明顯溫度拐點，且溫度拐點的深度隨時間的增加逐漸降低。當成孔12 h時孔內溫度已達到0℃(273.15 K)以下，此時可認為炮孔內部完全冰凍。在成孔2 h內深度在6m以內的炮孔溫度下降速率大，成孔2～24 h深度在2 m以內炮孔溫度變化劇烈。這是因為隨著時間的增加，外界地表溫度對炮孔內的溫度場影響逐漸減小，孔內滲流水因溫度的降低逐漸凍結為冰相，冰相凍結體的導熱系數大于空氣介質導熱系數，因此炮孔底部的低溫會迅速向孔外傳導，但由于炮孔頂部的地表溫度相對固定，所以造成在孔深6 m以內溫度顯著變化。

圖 7 炮孔成孔24 h內溫度隨深度變化圖Fig. 7 Temperature change with depth within 24 hours of blasthole formation

4.2 不同滲流水速度作用下炮孔內溫度分析

為進一步探討不同滲流水速度作用下炮孔內溫度的發展規律，在此模型上改變滲流水速度，分析流速為1 m/h、2 m/h、5 m/h、10 m/h時炮孔成孔3 h、8 h、12 h、24 h后孔內溫度變化，如圖8所示。

圖 8 不同滲流速度作用下溫度分布Fig. 8 Temperature distribution under different seepage velocities

由圖8可見：①孔內滲流水速度對炮孔內24 h后的整體溫度分布有較大影響。在一定流速范圍內，滲流水速度越快，炮孔內整體溫度越低。由于礦區底部溫度低且較為恒定，因此滲流水速度主要影響深度較淺的孔區溫度，對于深度較深的炮孔溫度影響不大。②隨著炮孔內溫度的降低，炮孔周邊巖體會隨滲流水的熱傳導作用溫度也相應降低。③從二維角度分析，炮孔成孔后不同速度的滲流水對礦區表面初始溫度的影響范圍≤2 m，對炮孔周邊巖體的影響范圍≤3 m，滲流速度越快，影響范圍越大。

分析不同滲流速度作用下孔內溫度隨時間變化規律及24 h時孔內溫度，不同滲流水速度作用下溫度對比如圖9所示，不同滲流水速度作用下24 h時孔內溫度如表2所示。

圖 9 不同滲流水速度作用下溫度對比Fig. 9 Temperature comparison under different seepage water velocities

表 2 不同滲流水速度作用下24 h時孔內溫度

由圖9可見：①在炮孔成孔后的任意時刻，滲流水速度越快，孔內溫度下降速度越快。②在不同滲流速度作用下，滲流速度越慢炮孔頂部溫度下降幅度越大，反之越小。當滲流速度1 m/h時，炮孔頂端溫度在24h內由277.92 K降為271.22 K，溫度下降幅度6.70K。滲流速度10 m/h時，炮孔頂端溫度在24 h內由267.16 K降為263.10 K，溫度下降幅度4.06 K。

由表2可見：①當滲流水為1 m/h時，孔內整體溫度都已達到0℃(273.15 K)以下，此溫度極有可能導致孔內裂隙水產生冰凍。②炮孔內最低溫度位于炮孔最深處且溫度恒定，當孔內裂隙滲流水速度高于5 m/h時24 h后炮孔內最大溫差將小于3 K，平均溫度為262.66 K，此時滲流水流速的增加對炮孔溫度變化影響不明顯。

4.3 不同初始地表溫度狀態下炮孔內溫度分析

由于不同地表溫度對炮孔深度較淺區域溫度影響較大，因此研究初始地表溫度對炮孔溫度的影響范圍及效果。依次改變地表初始溫度為0℃(273.15 K)及10℃(283.15 K)，對裂隙水滲流速度為3 m/h的場區溫度進行模擬，如圖10所示。對相同初始地表溫度狀態下不同時間炮孔內溫度變化進行對比，如圖11、圖12所示。

圖 10 不同初始地表溫度對炮孔溫度影響Fig. 10 The influence of different initial surface temperature on blasthole temperature

圖 11 地表初始溫度0℃(273.15 K)孔內溫度變化Fig. 11 The initial temperature of the ground surface is 0℃(273.15 K) and the temperature change in the hole

圖 12 地表初始溫度10℃(283.15 K)孔內溫度變化Fig. 12 The initial temperature of the ground surface is 10℃(283.15 K) and the temperature change in the hole

由圖10可見：①隨著初始地表溫度的上升，礦區深度在6 m內區域溫度明顯升高，深度6 m以上區域溫度變化不明顯，即初始地表溫度對礦區溫度影響限一定范圍。②從地表初始溫度為0℃(273.15 K)及10℃(283.15 K)溫度影響范圍對比可發現，地表溫度的升高會減弱滲流水對礦區溫度的影響范圍。

由圖11及圖12可見：①當礦區地表初始溫度為0℃(273.15 K)時，炮孔成孔3 h時整體溫度可達到零下，孔內平均溫度265.39 K，距離炮孔深度4 m以外時溫度變化相對于4 m以內較為明顯。炮孔成孔8 h后孔內整體溫度變化減緩，最高溫度267.39 K，孔內平均溫度264.44 K。12 h、24 h時炮孔平均溫度分別為264.14 K、263.74 K，溫度達到相對穩定值。當礦區地表初始溫度為10℃時，礦區地表溫度與礦區底部溫度梯度大，因此地表傳熱現象較為顯著。②從炮孔成孔后的4個時間節點可發現，隨著時間的增加，深度0～6 m處孔內溫度受外界環境影響大，孔深6 m以下區域溫度基本保持穩定。

5 現場試驗與模擬對比分析

試驗地點位于拉薩市墨竹工卡縣甲瑪露天礦海拔高度5120 m，天氣陰，溫度8～20℃，礦區頂平均氣溫6℃，最低溫度-2℃。

5.1 試驗儀器及使用

分布式光纖測溫系統(DTS)也稱為光纖測溫，主要原理以及喇曼(Raman)散射效應和光時域反射(OTDR)對光纖溫度的感度實現溫度監測。本試驗采用測溫儀器為SG-DTS-84U光纖測溫系統，測溫系統需使用熔纖機將分布式光纖與SG-DTS-84U光纖測溫系統串聯，待測前利用紅光筆檢測光纖通路，確認聯通后方可進行現場試驗。

5.2 試驗結果

通過炮孔成孔后24 h內現場試驗結果與模擬對比分析。SG-DTS-84U光纖測溫系統如圖13所示，炮孔內部圖如圖14所示，孔內滲流水溢出圖如圖15所示。

圖 13 SG-DTS-84U光纖測溫系統Fig. 13 SG-DTS-84U Fiber Optic Temperature Measurement System

圖 14 炮孔內部圖Fig. 14 Internal view of the blasthole

圖 15 孔內滲流水溢出Fig. 15 Seepage water overflow in the hole

通過現場試驗及COMSOL-Multiphysics有限元數值模擬，得到24 h內不同時間段炮孔內溫度變化情況，對比圖如圖16所示，孔內溫度誤差如表3所示。

圖 16 24 h內不同時間段炮孔內溫度對比圖Fig. 16 Comparison of the temperature in the blasthole at different time periods in 24 h

表 3 孔內溫度誤差

由圖16可見：①炮孔成孔后，最低溫度一般為炮孔底部，孔底低溫通過孔內滲流水逐漸向炮孔頂部擴散，在成孔2 h時0℃深度約為4 m，成孔8 h左右0℃深度約為2～3 m，成孔16 h左右0℃深度約為2 m，成孔24 h后孔內溫度可達到0℃(273.15 K)以下。②炮孔成孔后孔內溫度在8 h后下降速率明減慢并逐漸趨于穩定。

由表3可見：①實測溫度與模擬溫度在深度為0～4 m時溫差較大，最大溫差可達5 K，深度為4～18 m時誤差較小，一般不會高于2 K。這是因為在此深度的孔溫受外界環境影響較大，且實驗地點為高海拔地區，氣溫并不穩定，而炮孔內部溫度相對恒定，受外界因素干擾較小，因此接近地平面處孔溫會發生較大變化。

6 結論

以拉薩市墨竹工卡縣甲瑪礦區為研究對象，使用有限元模擬及現場試驗分析了高原高寒地區炮孔成孔后孔內溫度場變化規律，所得結論如下：

1)炮孔成孔后溫度變化主要受地層溫度影響，導致孔內積水凍結，且炮孔溫度隨深度整體呈下降趨勢，初始地表溫度對炮孔上部的影響趨于一定范圍(≤6 m)，與炮孔上部(0～5 m)區域在受空氣對流熱通量耦合作用影響致使孔內積水未達到凍結臨界狀態。

2)孔內24 h后溫度變化主要由滲流水速度決定，但當流速達到溫度影響流速閾值(5 m/h)時，滲流水速度對炮孔內溫度變化影不再顯著。

3)炮孔成孔8 h后溫度場會達到相對穩定值，裝藥更為合理，若因特殊原因無法在此時間內裝藥，筆者會在今后的研究方案中提出應對措施。

(感謝成遠礦業開發股份有限公司及羅乃鑫、姚毅、唐玲彪、白宇在科研現場試驗期間給予的支持！)