張 澤， 王述紅， 楊天嬌， 張雨濃

(東北大學 資源與土木工程學院， 遼寧 沈陽 110819)

我國東北部以及西部這些高緯度或高海拔地區(qū)有一半以上的國土屬于寒區(qū)[1].寒區(qū)隧道在修建過程中由于溫度降低圍巖會受到凍脹的影響，除了巖土體孔(裂)隙中的原位水發(fā)生凍結(jié)外，巖土體中的水分遷移使得水分不斷向凍結(jié)鋒面遷移，促使凍脹更加明顯.因此，隧道巖體中的凍脹問題歸根結(jié)底為土體中的溫度水分遷移問題.凍融狀態(tài)轉(zhuǎn)變過程中水分遷移導致土體中水分占比變化引起土體熱參數(shù)發(fā)生改變；融化界面的冰晶會對周圍的水分產(chǎn)生吸引作用從而導致界面水頭的產(chǎn)生；由于水冰相變土體中水分含量降低而冰含量相應升高因此阻礙水的擴散，這些會對凍結(jié)過程造成深遠影響；同時相變釋放的熱量也會對凍結(jié)產(chǎn)生較大影響；即使在凍結(jié)溫度時土體中仍然存在一部分未凍水，且其含量與負溫保持動態(tài)平衡關系.

隧道凍害過程是溫度-滲流耦合問題，長期以來，各國學者建立了不同的溫度-滲流耦合模型對寒區(qū)隧道的凍融進行研究[2-4].Liu等[5]采用熱力學、流體力學、彈塑性力學基本理論，建立了非飽和多孔介質(zhì)材料低溫THM耦合模型，對凍脹量的分布、凍脹深度以及應力場、溫度場的分布進行模擬計算分析.譚賢君等[6]通過數(shù)值模擬對提出的“三區(qū)域”理論進行分析，但是對于水分遷移以及水冰相變的阻礙作用沒有做出解釋.Yamabe等[7-8]對飽和的砂巖和凝灰?guī)r進行了實驗并對低溫水冰相變下的飽和巖石多次耦合過程進行了研究.徐光苗等[9]首次對低溫巖石進行研究，建立了低溫巖石的THM三場耦合方程，并采用有限元程序進行模擬.

本文基于混合物理論，考慮水分遷移以及水冰相變的影響，以及土體孔隙冰對未凍水的阻滯作用對原有模型進行改進并應用于非飽和土的低溫水熱耦合問題中，并使用COMSOL軟件[10]進行模塊開發(fā)與文獻[11]中經(jīng)典實驗進行數(shù)值驗證，證明本文所建立的溫度-滲流耦合數(shù)學模型的適用性.然后以西藏自治區(qū)米林隧道為工程實例，對考慮水分遷移與不考慮水分遷移的溫度場、水分場進行了模擬分析，研究成果能較真實地反映寒區(qū)富水隧道凍害現(xiàn)象發(fā)生過程，具有一定的參考價值.

1 寒區(qū)隧道水熱耦合

為了建立更符合寒區(qū)隧道凍脹問題的溫度-滲流耦合模型，在現(xiàn)有經(jīng)典Harlar水熱模型基礎上，將水分遷移和水冰相變考慮其中建立更為全面的巖土介質(zhì)凍結(jié)過程的溫度-滲流耦合模型，為寒區(qū)隧道的凍害問題提供參考.

1.1 基本假設

1) 土介質(zhì)為均質(zhì)各向同性孔隙介質(zhì)，由未凍水、冰、巖體骨架組成，凍結(jié)過程中不發(fā)生變形.

2) 巖土介質(zhì)中水分遷移，忽略空氣對水分遷移的貢獻.

3) 只考慮水冰相變的潛熱和熱傳導過程.

4) 不考慮由于冰透鏡體等冰體形成造成的上部已凍土克服下部水分重力對水分的抽吸作用.

1.2 溫度場控制方程

根據(jù)傅里葉定律和能量守恒定律，在溫度場控制方程中的比熱容、導熱系數(shù)都為含水率的函數(shù)，結(jié)合凍土的相變條件，單元體內(nèi)的體積熱容可描述為單元體的冰水相變所產(chǎn)生或吸收的相變熱，由此可得

(1)

從而可得溫度場的控制方程：

(2)

式中：Ca為等效體積熱容量，kJ/(m3·℃)；t為土體瞬態(tài)溫度,℃；τ為時間,s；λ為熱傳導系數(shù)[W/(m·℃)]；ρw,ρi分別為水、冰密度,kg/m3；θi為冰體積含量；L為相變釋放潛熱，取值334 kJ/kg；為哈密頓算子.

1.3 水分場控制方程

在寒區(qū)隧道工程圍巖中即使在凍結(jié)點溫度以下仍然存在沒有凍結(jié)的水分，其遵循達西定律進行遷移：在單元體內(nèi)考慮孔隙冰對未凍水的阻滯作用，非飽和土中的未凍水遷移微分方程即為水分場耦合溫度控制方程：

(3)

式中,D(θu)，K(θu)分別為非飽和土擴散系數(shù)和導水率，可見其都是隨著未凍水含量發(fā)生變化.

凍土中水的擴散率D(θ)計算：

(4)

式中：k(θ)為非飽和土的滲透率,m/s；c(θ)為比水容量,1/m，可由滯水模型確定.

在凍土中，孔隙冰含量和孔隙水含量存在動態(tài)變化，冰的增加必占據(jù)孔隙從而減少水分遷移的通道.為了體現(xiàn)與常溫下非飽和土水分遷移的區(qū)別，考慮了負溫下孔隙冰的阻抗作用，引入阻抗系數(shù)的概念.I為阻抗系數(shù).

I=10-10θi，

(5)

k(θ)=ksSl(1-(1-S1/m)m)2,

(6)

c(θ)=a0m/(1-m)·S1/m·(1-S1/m)m.

(7)

其中：l為本構參數(shù)；S為凍土相對飽和度，定義為

(8)

式中,θr和θs分別代表土的剩余含水量和飽和含水量.

2 關鍵耦合參數(shù)描述

1) 等效熱容:

(9)

2) 等效熱傳導系數(shù).導熱系數(shù)的取值有多種加權方法，本文采用指數(shù)加權方法，可表示為

λ(θ)=(λs)θs(λl)θl(λi)θi.

(10)

3) 未凍水含量.為真正建立水分場與溫度場的耦合方程，需要引入θ,θi與t之間的聯(lián)系方程.徐學祖等[12]基于試驗數(shù)據(jù)，給出了凍土中未凍水含量與溫度的經(jīng)驗關系式：

(11)

式中：tf為土體凍結(jié)溫度,℃；w0為土體初始含水量,%；wu為負溫度為t時的未凍水質(zhì)量分數(shù),%；B是與土類和含鹽量確定的常量.按照徐學祖等[12]凍土中水分遷移的實驗研究中的測定方法，砂土、粉土、黏土中的B值為0.61，0.47，0.56.

4) 固液比.凍土中孔隙冰體積與未凍水體積之比，記為Bi:

(12)

因此，凍土中孔隙冰、未凍水和溫度的聯(lián)系方程為

θi=Bi(t)·θu.

(13)

3 水熱力模型的試驗驗證

文獻[11]中的一維土柱凍結(jié)試驗可以對所建立的模型正確性和可行性進行驗證，該試驗為封閉系統(tǒng),如圖1所示.

土樣由非飽和土組成，無外載凍結(jié)方式為從上至下，試驗參數(shù)見表1.

表1 試驗參數(shù)

文獻[11]中給出了試驗的初始溫度見圖2，土柱上邊界和下邊界溫度的變化值見表2.

表2 試驗上下邊界溫度變化

本文計算的物理參數(shù)如表3所示.

表3 計算物理參數(shù)

根據(jù)文獻中的實驗結(jié)果分別對528 min和2 830 min 的溫度變化進行計算，與文獻[11]中的試驗數(shù)值進行了對比，如圖3所示，驗證結(jié)果基本吻合，證明了模型的正確性.

4 工程實例

4.1 工程概況及計算參數(shù)

米林隧道址山勢雄偉，區(qū)內(nèi)高點位于隧道軸線左側(cè)山脈，海拔高度為4 230 m，最低點位于隧道出口寬谷地帶，標高為2 940 m，典型的高海拔隧道，最大埋深為1 200 m.項目所在地區(qū)歷年年平均氣溫4.3 ℃，歷年極端最高氣溫為21 ℃，歷年極端最低氣溫為-25 ℃，歷年最冷月平均氣溫為-7.6 ℃，按對隧道工程影響的氣候分區(qū)，屬寒冷地區(qū).米林隧道為高海拔富水寒區(qū)大長隧道.

根據(jù)地質(zhì)勘查資料，選取隧道的上覆山體和隧道下一部分山體進行計算，取其硐深10 m處的斷面為研究對象，該斷面處的地質(zhì)條件為：地質(zhì)成分主要為塊石土、粉土.根據(jù)隧道的埋深以及斷面尺寸，建立隧道均勻孔隙介質(zhì)簡化計算模型，如圖4所示.

1) 米林隧道溫度場初始條件和邊界條件.

林芝縣近五年月平均氣溫的變化情況：林芝縣全年的平均氣溫為2.3 ℃，最冷月為一月,平均氣溫為-9.8 ℃，最熱月為7月,平均氣溫為11.8 ℃；極端最高氣溫為27 ℃，極端最低氣溫為-13 ℃，由于海拔高，溫度垂直變化明顯，海拔每升高100 m溫度下降0.6 ℃.

將五年的每月平均溫度繪制如圖5所示,大致認為一年內(nèi)的溫度變化規(guī)律大致呈現(xiàn)周期性變化，可寫成正弦或余弦函數(shù)的形式:

(14)

實際中隧道頂部和隧道內(nèi)部的大氣溫度變化會引起圍巖中的溫度-滲流變化.因此選取隧道斷面施工期的地面(ab邊)及內(nèi)部(efg邊)溫度為大氣溫度變化，見式(14).根據(jù)對稱性ae,gd邊為絕熱邊界，同樣bc,cd邊也為絕熱邊界.由于為高海拔隧道，需要考慮高度對于溫度的影響，隧道頂部溫度要低于隧道圍巖溫度如圖6所示.根據(jù)施工現(xiàn)場提供的隧道施工階段溫度值，隧道圍巖初始溫度為-0.98 ℃.

本文主要研究低溫狀態(tài)下隧道中水熱耦合效應，因此選取本年11月到次年4月這六個月進行研究，溫度取值變化趨勢如圖6所示.

2) 滲流場初始邊界條件與邊界條件.

bc邊為水頭壓力邊界，若將d點作為原點進行計算，水壓力在豎直方向上呈線性變化，總水頭等于基質(zhì)吸力水頭與位置水頭之和，因此取壓力水頭為式(15)所示，初始含水量根據(jù)現(xiàn)場提供數(shù)據(jù)取值為12%.

(15)

式中：p為水頭壓力，Pa；h為高度，m.

根據(jù)對稱性ae，gd為不透水邊界，襯砌與圍巖的接觸邊efg壓力水頭取值為0，cd邊為固定邊界，壓力為0，忽略襯砌的滲透性，初始水壓力取值為0.

由地質(zhì)條件可知，米林隧道計算所取斷面主要為粉土，并噴射60 cm厚的C25混凝土襯砌.主要計算物理參數(shù)如表4所示，其他參數(shù)如表3所示.

表4 模型計算主要參數(shù)

4.2 米林隧道水熱耦合模型分析

1) 溫度場時空分布.林芝地區(qū)在11月份開始進入零下溫度，并持續(xù)降溫.圖7顯示了隧道整體由11月份到次年4月份經(jīng)歷的溫度變化狀況.根據(jù)當?shù)氐臍鉁刈兓?guī)律，隧道邊界氣溫持續(xù)降低，隧道頂部邊界氣溫由-0.82 ℃降低到-9 ℃，隧道內(nèi)部邊界溫度由-0.74 ℃降低到-11.11 ℃，導致隧道整體溫度的下降，下降趨勢為由頂部向底部，由洞口向內(nèi)部延伸擴展.

距離邊界越近的巖體溫度下降越快，這是因為邊界帶范圍內(nèi)的溫度梯度大，在導熱系數(shù)相同的情況下，熱量傳遞量大；此外再考慮凍結(jié)狀態(tài)下水冰相變的情況， 水相變?yōu)楸鶗r比熱容會降低而導熱系數(shù)會增大，由此也增加熱量的傳輸速率.在3月和4月溫度有明顯的回升.

2) 未凍水含量時空分布.由上述溫度時空分布的分析知，在11月到次年2月的120 d中隧道頂部及隧道內(nèi)部經(jīng)歷降溫的過程，隧道巖體存在局部降溫，降溫范圍可以從圖8a～圖8d觀察得出，紅色區(qū)域持續(xù)減小，說明此區(qū)域存在持續(xù)降溫狀態(tài).

由圖8e,圖8f分析可得，由于隧道頂部及內(nèi)部溫度升高導致低溫時已經(jīng)凍結(jié)成冰的固態(tài)水轉(zhuǎn)化為液態(tài)水， 顯示了隧道頂部邊界和隧道邊界由模擬三月時逐漸轉(zhuǎn)變?yōu)檎郎?，說明邊界進入了融化狀態(tài)，導致隧道頂部邊界和隧道內(nèi)部邊界未凍水含量升高.在隧道頂部和隧道內(nèi)部未凍水含量由120 d時的5.3%迅速增加到180 d時的9.8%，這對隧道的穩(wěn)定性是十分不利的.

3) 冰含量時空分布.由圖9a～圖9d可以看出，隧道巖體的固態(tài)冰含量持續(xù)升高，這是因為隧道巖體持續(xù)降溫，液態(tài)水被凍結(jié)為固態(tài)水，從而升高了固態(tài)冰的含量，由上述溫度時空分布的分析知，在11月至2月隧道頂部以及隧道內(nèi)部經(jīng)歷降溫過程，隧道巖體中存在局部降溫，降溫的范圍可以從圖9a～圖9d觀察出來，藍色區(qū)域(固態(tài)冰含量5%～10%)持續(xù)減小，由藍綠色區(qū)域的擴張代替(固態(tài)冰含量10%～25%)，說明此區(qū)域經(jīng)歷著持續(xù)降溫狀態(tài)，并且含冰量峰值出現(xiàn)在1月份.

由圖9e，圖9f分析可得，由于隧道頂部以及隧道內(nèi)部的升溫導致已經(jīng)凍結(jié)成冰的固態(tài)水轉(zhuǎn)為液態(tài)水，因此隧道在3月份開始內(nèi)部的凍結(jié)冰開始融化，在4月份時隨著溫度的持續(xù)升高，其含冰量也隨之降低，由凍結(jié)時的25%降低到15%，可以推測隨著溫度的持續(xù)升高，隧道巖體中的冰含量將一直減小直到全部融化為液態(tài)水.巖體中的水冰相變和冰體的凍結(jié)融化對于隧道整體的穩(wěn)定性產(chǎn)生重要的影響.

4) 水分遷移對溫度場的影響分析.為與上述已經(jīng)考慮水分遷移的模型對比，應用式(16)建立不考慮水分遷移的數(shù)值模型，只考慮隧道巖體整體的熱遷移，不考慮冰水相變.

(16)

通過對比圖10和圖11可見，未考慮水分遷移的溫度場中熱傳導速度較快，圖10的紅色區(qū)域較圖11的紅色區(qū)域減小的速率大，其黃綠色區(qū)域(-4 ℃～-7 ℃)擴張速度大，這是因為考慮水分遷移的熱傳導過程受到的水相變成冰釋放潛熱的影響，相變潛熱提供的熱量阻止了低溫的傳遞，而且水的比熱容一般是巖體的比熱容的4倍左右，說明巖石升高同等溫度要消耗4倍的熱量，也在一定程度上阻礙了隧道巖體中的溫度升高過程.

由3月份溫度場對比分析可以看出，考慮水分遷移的隧道頂部以及隧道內(nèi)部溫度場與不考慮水分遷移的溫度場相比，前者升溫過程也比較緩慢，分析原因，導致這種情況主要有兩點：一是由于水的比熱容較大，升高相同溫度消耗的熱量較巖體的多；二是隨著溫度周期性變化，由于巖體升溫導致巖體內(nèi)的冰含量減少，而冰融化成水是吸收熱量的過程，要消耗熱能，相對的水成冰要放熱，冰成水要吸熱，所以考慮水分遷移的隧道巖體在升溫過程中的傳熱較慢.

5 結(jié) 論

1) 在原有水熱耦合模型基礎上，推導出考慮水冰相變和水分遷移的水熱耦合方程，從而建立了更全面更符合實際的水熱耦合模型.利用COMSOL軟件開發(fā)新的模塊與一維土柱凍結(jié)實驗進行反演驗證，數(shù)值模擬得到的溫度場與實驗結(jié)果基本吻合較好，進而驗證了模型的正確性.

2) 利用水熱耦合模型對米林隧道11月到次年4月的溫度場和水分場進行模擬分析.結(jié)果表明：隧道頂部邊界氣溫由-0.82 ℃降低到-9 ℃，隧道內(nèi)部邊界溫度由-0.74 ℃降低到-11.11 ℃，并在3月份隧道溫度回升；同時隧道中冰含量在1月份達到峰值，在3月份開始減小，表明在隧道中存在明顯的凍融現(xiàn)象，極易產(chǎn)生凍害問題.

3) 對比分析水分遷移對溫度場的影響，考慮水分遷移的隧道頂部以及隧道內(nèi)部溫度場升溫過程比較緩慢，證明固態(tài)冰的形成或融解時所產(chǎn)生的相變潛熱對隧道中溫度場的分布影響遠遠大于液態(tài)水靠重力遷移造成的熱對流傳熱.