寒區(qū)隧道圍巖水熱耦合數(shù)值分析

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

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

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

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

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

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

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

1.1 基本假設(shè)

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

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

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

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

1.2 溫度場(chǎng)控制方程

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

(1)

從而可得溫度場(chǎng)的控制方程：

(2)

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

1.3 水分場(chǎng)控制方程

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

(3)

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

凍土中水的擴(kuò)散率D(θ)計(jì)算：

(4)

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

在凍土中，孔隙冰含量和孔隙水含量存在動(dòng)態(tài)變化，冰的增加必占據(jù)孔隙從而減少水分遷移的通道.為了體現(xiàn)與常溫下非飽和土水分遷移的區(qū)別，考慮了負(fù)溫下孔隙冰的阻抗作用，引入阻抗系數(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為本構(gòu)參數(shù)；S為凍土相對(duì)飽和度，定義為

(8)

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

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

1) 等效熱容:

(9)

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

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

(10)

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

(11)

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

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

(12)

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

θi=Bi(t)·θu.

(13)

3 水熱力模型的試驗(yàn)驗(yàn)證

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

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

表1 試驗(yàn)參數(shù)

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

表2 試驗(yàn)上下邊界溫度變化

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

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

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

4 工程實(shí)例

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

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

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

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

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

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

(14)

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

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

2) 滲流場(chǎng)初始邊界條件與邊界條件.

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

(15)

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

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

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

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

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

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

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

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

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

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

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

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

(16)

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

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

5 結(jié) 論

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

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

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