多圈管凍結壁水熱時空分布數值模擬研究

羅陽， 王浩， 倪斌

（1.安徽理工大學土木建筑學院，安徽 淮南 232001；2.安徽壽縣農村商業銀行股份有限公司，安徽 淮南 232281）

0 引言

人工凍結法是一種使用人工制冷技術對地層進行人工制冷，使得天然巖土變成凍土，形成人工凍結壁，由于凍土的高抗壓抗剪性能以及良好的阻水性能，在人工凍結壁內側形成了良好的施工環境，有利于地下工程的施工［1］。但人工凍結壁［2］的形成過程是一個復雜的多物理場過程，其中主要包括溫度場、水分場、應力場的相互作用過程，這個復雜的多物理場過程存在于同一圈不同凍結孔之間，也存在于不同凍結圈之間以及凍結壁與天然巖土層之間。

目前國內人工凍結法的應用主要在一些大型礦井的開挖以及隧道掘進等工程上，人工凍結法根據主凍結圈+輔助凍結圈的凍結形式，從最初的單圈管凍結發展到雙圈管、三圈管、四圈管。國內學者［3］對多圈管凍結的研究主要為室內模型試驗，通過模型試驗研究多圈管凍結壁的形成過程，還有學者通過數值仿真軟件研究了多圈管凍結壁溫度場發展規律。汪仁和等［4］開展了大型的人工多圈管凍結模型試驗，研究了多圈管凍結壁的形成和融化過程中的溫度場、水分場、應力場，提出了導水系數是溫度梯度的函數的結論。李棟偉等［5］利用ADINA 軟件進行了多圈管凍結瞬態溫度場的有限元數值仿真，獲得了凍結壁的有效厚度、凍結壁的平均溫度，與現場數據較為吻合。認為有限元溫度場數值仿真可以較好地預測多圈管凍結的溫度場。陳軍浩等［6］進行了多圈管凍結壁的形成與融化研究，進行了多圈管模型優化設計的室內模型試驗。

文中通過對COMSOL Mutiphysics 軟件進行了二次開發，實現了多圈管凍結壁的水熱耦合數值仿真，通過與現場數據對比，表明正凍非飽和土水熱耦合模型可以用于描述多圈管凍結過程中的水熱分布。通過對數值模擬結果進行分析，進行了多圈管凍結壁各特征點溫度、冰水質量、凍結時間關系分析，以及主界面水熱時空分布分析。

1 正凍非飽和土水熱耦合模型

1.1 熱量遷移方程

基于能量守恒定律和熱傳導理論，將冰水相變作為內熱源處理，建立如下熱量遷移方程［7］：

式中，ρ 和ρi為土的密度和冰的密度，kg/m3；C（θ）為體積熱容，J/（kg℃）；λ（θ）為導熱系數，W/（m℃）；T為土體的瞬態溫度，℃；t為時間，s；△為微分算子；θ為體積含水率；θi為孔隙冰體積含量；L為相變潛熱，通常取值為334.5kJ/kg。

1.2 水分遷移方程

根據Richard方程［8］，基于質量守恒定律和達西滲流理論，非飽和凍土中的未凍水的水分遷移方程：

式中，D，K 分別為土擴散系數和導水率，均為未凍水含量θu的函數。

1.3 聯系方程

通過對徐學組建立的未凍水預測模型變形，獲得固液比與溫度的關系［9］。

冰的體積含量θi三者的聯系方程：

以上，熱量遷移方程、水分遷移方程、聯系方程3個方程構成的方程組建立了正凍非飽和土的水熱耦合模型。

2 多圈管凍結水熱耦合模擬

2.1 幾何建立

將水熱耦合數值模擬理論應用于皖北地區某礦副井的多圈管凍結方案，幾何模型為副井平面模型的1/4，凍結孔、凍結孔圈尺寸、和凍結孔位置均按實際工程取值。

2.2 邊界條件及物理參數

土體所有外邊界假定為絕緣邊界，內部凍結孔溫度場邊界取實際現場循環液去、回路溫度的平均值，內部凍結孔溫度邊界曲線如圖1所示。水分場邊界假定為絕緣邊界，凍結壁周圍初始含水率為28%。深厚黏土熱量特征參數、土水分特征參數見表1、表2。

圖1 凍結孔溫度邊界

表1 深厚黏土熱量特征參數

表2 深厚黏土土水分特征參數

2.3 數值模擬結果及分析

2.3.1 凍結壁厚（180d）。

當達到180d，內邊界與輔助凍結圈間距大于凍結壁外邊界與主凍結圈的間距，這是因為凍結壁內側土體在主孔包裹下始終封閉狀態，外側則與外圈的土體發生熱交換。

圖2 不同時刻凍結壁發展

2.3.2 溫度場

圖3演示了凍結孔附近土體溫度在不同時間內的發展情況：凍結孔降溫到指定溫度（30d）→主凍結孔圈降溫交圈（60d）→輔助孔圈降溫交圈（90d）→主、輔助孔溫度分布呈規則圓環（120d）→溫度分布均勻且有所回升（180d）。

圖3 不同凍結時刻溫度場分布

在整個過程中可以看出，整個過程中中心土體的溫度始終降低直到負溫。在凍結的第180d時，主凍結圈和輔助凍結圈間的土體溫度有所上升，溫度在-15℃左右，這與循環液溫度上升，凍結進入養護階段有關。

2.3.3 水分場

（1） 未凍水含量。圖4 演示了凍結孔周圍土體含水量隨著時間的變化過程，在整個過程可以很明顯的看到：①主、輔助孔中間土體含水量最先降低，從開始的14%（30d）降低至8%（60d）直到降至最低含量6%（90d、鋸齒狀），在凍結的第180d，主、輔助孔之間未凍水含量低于6%土體范圍逐漸減小，呈魚刺狀；②主、輔助孔邊界土體未凍水含量低于10%的邊界土體呈規則環形向內和向外發展。

圖4 不同時刻未凍水分布

（2） 含冰量。根據圖5 不同時刻冰水分布圖分析可知：①在凍結的第30d，凍結孔附近土體含冰量達到24%左右，僅在凍結孔附近很小的范圍內有土體冰水含量達到30%以上，在主凍結圈外圈有一環帶含冰量小于18%，主凍結圈和輔助凍結圈間的土體含冰量低于10%，呈魚刺狀分布；②在凍結的第60d，主凍結圈和輔助凍結圈間低于18%的土體魚刺狀分布消失，在主凍結圈外和輔助凍結圈內土體呈冰水含量低于18%的環形帶，隨著凍結的持續，低于18%的環形帶向內和向外移動，在主凍結圈外側形成一條冰水含量為25%的均值環形帶；③整個凍結過程中，含冰量大于30%的土體僅存在在凍結孔附近的較小范圍內。

圖5 不同時刻冰水分布

3 數值結果與現場數據對比

3.1 溫度變化

從圖6可以看出：①1#、3#測點為主凍結圈外1m處測點，數值計算結果為二者的平均值，變化規律一致，2#測點為主、輔助凍結圈中間測點，數值模擬結果與實測結果在凍結的40～60d 完全重合，曲線變化規律一致；②在試驗后期溫度均呈上升趨勢，數值計算結果低于實測結果，這和計算模型中忽略土體自身溫度隨季節性變化相關。

圖6 1#、3#、2#溫度變化

基于數值模擬的計算結果與現場實際情況較吻合，數值模擬的計算過程能夠清楚描述多圈管凍結過程中的熱時空分布。

3.2 水-熱-時間關系

根據圖7分析可得：①主、輔助凍結圈中間土體與輔助凍結孔中間土體溫度下降過程中出現明顯的相變潛熱釋放過程；②凍結180d 后，土體溫度由低到高排序：主凍結圈與輔助凍結圈位置、主凍結孔中間位置、輔助凍結孔中間位置、井幫位置、主凍結圈外1m處。前三的溫度值相差不大，均為-15℃左右，井幫位置溫度為-10℃左右，主凍結圈外1m 處溫度為-6℃左右；③各特征點未凍水含量隨時間變化關系為：主、輔助凍結圈位置、主凍結孔中間位置、輔助凍結孔中間位置三處未凍水較早時間內達到平衡狀態主凍結圈外1m 處未凍水全過程處于緩慢下降階段，井幫位置處未凍水變化經歷慢-快-慢三階段；④所有特征點的冰水含量變化均為先減小后增大。開始增大時該處土體溫度為凍結溫度，即土體降至凍結溫度時，土體冰水含量開始增大。主凍結圈外1m處土體冰水質量含量最大，達到23.5%左右，主凍結中間位置土體冰水質量含量最小，為19.5%左右，其他三處為20%左右。

圖7 不同特征點水-熱-時間關系

3.3 主、界面溫度、冰水含量分布

從圖8可以看出：①凍結初期，從井筒中心到凍結壁外側，主面溫度：高、低、高、低、高，隨著凍結的持續，在主凍結圈和輔助凍結圈之間的土體溫度逐漸變得平緩，最終達到和凍結孔一致的溫度；②界面溫度分布與主面分布變化一致，但界面溫度較主面溫度分布平緩。

圖8 主、界面溫度分布

從圖9可以看出：①主面和界面的水分分布為高、低、高、低、高、低、高。隨著凍結的持續，在凍結壁的外側出現冰水含量二次高點，這是由于凍結壁向外發展較緩慢，在凍結壁邊緣出現水分集聚現象；②同半徑距離，主面的冰水含量高于界面，冰水含量大于28%的土體范圍僅存在各凍結孔附近直徑約0.3m 范圍內；③凍結壁兩側的低含水帶隨凍結壁的增大而向兩側移動；④主、輔凍結圈中間的土體水分在凍結初期呈減小狀態，后隨著凍結的持續，水分有所增大，在凍結60d后無明顯變化。

圖9 主、界面冰水分布

4 結語

（1） 基于正凍非飽和土水熱耦合模型，進行COMSOL Multiphysics 軟件二次開發用于水熱耦合計算，能夠較好地描述多圈管凍結過程中的水熱時空分布，對實際工程有一定指導意義。

（2） 主凍結圈+輔助凍結圈形式的雙圈管凍結，主凍結圈首先發生交圈，主凍結圈完全交圈后，輔助凍結孔先于主凍結圈交圈再與臨近輔助凍結孔交圈，在交圈的過程中，主、輔凍結圈間會出現局部未凍倉。

（3） 未凍水含量較低的土體首先出現在主凍結孔之間，隨著凍結的持續，在主、輔凍結圈間呈鋸齒狀分布，后由于凍結從積極期進入維護期，鋸齒狀分布變為魚刺狀分布，范圍減小。

（4） 凍結壁范圍內各特征點處的冰水含量均經歷先減小后增大兩個階段，轉折點出現在該處土體溫度處于凍結溫度，即土體凍結后，土體水分含量開始增大。

（5） 主面溫度分布較界面溫度分布平緩，在凍結維護期，無論主面還是界面，主、輔凍結圈間的土體溫度均與凍結孔處溫度一致。

（6） 同半徑范圍，主面的冰水含量要高于界面，隨著凍結的持續，在凍結壁外側出現冰水含量次高峰，凍結壁兩側的低含水帶隨凍結壁的增大而向兩側移動，主、輔凍結圈中間的土體水分在凍結初期減小，后隨著凍結的持續，水分有所增大，在凍結60d后無明顯變化。