羅陽， 王浩， 倪斌

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

0 引言

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

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

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

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

1.1 熱量遷移方程

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

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

1.2 水分遷移方程

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

式中，D，K 分別為土擴散系數(shù)和導(dǎo)水率，均為未凍水含量θu的函數(shù)。

1.3 聯(lián)系方程

通過對徐學(xué)組建立的未凍水預(yù)測模型變形，獲得固液比與溫度的關(guān)系［9］。

冰的體積含量θi三者的聯(lián)系方程：

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

2 多圈管凍結(jié)水熱耦合模擬

2.1 幾何建立

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

2.2 邊界條件及物理參數(shù)

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

圖1 凍結(jié)孔溫度邊界

表1 深厚黏土熱量特征參數(shù)

表2 深厚黏土土水分特征參數(shù)

2.3 數(shù)值模擬結(jié)果及分析

2.3.1 凍結(jié)壁厚（180d）。

當(dāng)達(dá)到180d，內(nèi)邊界與輔助凍結(jié)圈間距大于凍結(jié)壁外邊界與主凍結(jié)圈的間距，這是因為凍結(jié)壁內(nèi)側(cè)土體在主孔包裹下始終封閉狀態(tài)，外側(cè)則與外圈的土體發(fā)生熱交換。

圖2 不同時刻凍結(jié)壁發(fā)展

2.3.2 溫度場

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

圖3 不同凍結(jié)時刻溫度場分布

在整個過程中可以看出，整個過程中中心土體的溫度始終降低直到負(fù)溫。在凍結(jié)的第180d時，主凍結(jié)圈和輔助凍結(jié)圈間的土體溫度有所上升，溫度在-15℃左右，這與循環(huán)液溫度上升，凍結(jié)進(jìn)入養(yǎng)護(hù)階段有關(guān)。

2.3.3 水分場

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

圖4 不同時刻未凍水分布

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

圖5 不同時刻冰水分布

3 數(shù)值結(jié)果與現(xiàn)場數(shù)據(jù)對比

3.1 溫度變化

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

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

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

3.2 水-熱-時間關(guān)系

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

圖7 不同特征點水-熱-時間關(guān)系

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

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

圖8 主、界面溫度分布

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

圖9 主、界面冰水分布

4 結(jié)語

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

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

（3） 未凍水含量較低的土體首先出現(xiàn)在主凍結(jié)孔之間，隨著凍結(jié)的持續(xù)，在主、輔凍結(jié)圈間呈鋸齒狀分布，后由于凍結(jié)從積極期進(jìn)入維護(hù)期，鋸齒狀分布變?yōu)轸~刺狀分布，范圍減小。

（4） 凍結(jié)壁范圍內(nèi)各特征點處的冰水含量均經(jīng)歷先減小后增大兩個階段，轉(zhuǎn)折點出現(xiàn)在該處土體溫度處于凍結(jié)溫度，即土體凍結(jié)后，土體水分含量開始增大。

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

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