鮮少華，李清鵬，羅振先

(1.武漢市政工程設計研究院有限責任公司， 湖北 武漢 430023；2.河南大學 土木建筑學院， 河南 開封 475004)

當季節(jié)凍土區(qū)地下水埋深較淺時，在冬季凍結過程中水分會在溫度梯度的作用下遷移至路基土體中，引起凍脹變形，繼而產生路面開裂破壞；融化過程中，大量的遷移水不易消散，造成土體強度降低，進而引起路面沉陷等病害。溫度梯度的變化導致水分遷移的過程即是一種水熱耦合的過程。水熱耦合理論一直是國內外學者的研究熱點，Harlan[1]通過對比分析非飽和土與部分凍結土中的水分遷移規(guī)律，最早建立了描述水熱遷移的數(shù)學模型。隨后，Taylor等[2]、Jame等[3]、Newman等[4]不同程度地改進和發(fā)展了水熱遷移模型。Mu等[5]引入固體力學理論，在計算體應變時考慮了水分凍脹引起的土體應變，提出了水熱力三場耦合模型。毛雪松等[6]利用水熱力耦合模型計算了青藏公路唐南段路基溫度場、濕度場和應力分布規(guī)律，得到的結果能夠體現(xiàn)道路凍脹的主要誘因。張玉芝等[7]聯(lián)合土體凍結過程中的熱傳導方程和熱彈性力學理論建立了高鐵路基的變形和內力計算模型，預測了路基服役期內在凍融作用下可能產生的凍脹量和融沉量。邰博文等[8]在凍土水熱耦合方程的基礎上通過土體凍脹率與含冰量的關系改進了土體凍脹模型，并利用寒區(qū)路基變形監(jiān)測數(shù)據(jù)驗證了模型的準確性。Zheng等[9]將一維凍脹模型擴展到了二維情形并考慮了凍結速率和約束應力的影響。Li等[10-11]在水熱力三場耦合模型將凍融循環(huán)作用對土體的影響考慮了進去，并利用室內試驗和現(xiàn)場監(jiān)測進行了驗證。褚志成等[12]借助溫度場-應力場耦合模型研究了多年凍土區(qū)邊坡在凍融作用下的局部穩(wěn)定性。呂龍等[13]在Harlan水熱耦合模型的基礎上考慮水分遷移和冰水相變，認為水分遷移是導致基坑工程冬季凍脹破壞的主要原因之一。

綜上，國內外諸多研究者在土體凍脹機理和水熱耦合理論方面已經進行了大量研究工作，并且取得了一定的成果，但上述研究均未考慮有外界水分持續(xù)補給的情況。鑒于此，基于現(xiàn)有的水熱耦合理論，針對季節(jié)凍土區(qū)路基填土的水熱遷移問題，采用室內試驗和數(shù)值計算的手段，研究有外界持續(xù)補水條件下土體內部的溫度場和濕度場變化規(guī)律。

1 水-熱耦合理論

1.1 基本假定

基于現(xiàn)有理論[14]，作以下基本假定：(1) 土體均勻連續(xù)，屬于各向同性體；(2) 土中無鹽分影響；(3) 凍結區(qū)和未凍區(qū)水分僅以液態(tài)形式進行遷移；(4) 土顆粒不可壓縮；(5) 已凍土和未凍土為彈性體；(6) 無溶質遷移；(7) 水分遷移符合Darcy定律。

1.2 溫度場方程

忽略對流傳熱和熱輻射，考慮土體中冰水相變的熱傳導方程表示為：

(1)

1.3 水分場基本方程

利用水動力學模型來描述土體凍結過程中水分的遷移過程。對于飽和土和非飽和土，水分遷移方程均可表示如下：

(2)

式中：θw為水分總體積含量，即θw=θu+θiρi/ρw，其中θu為未凍水的體積含量；ρw為水的密度；kx，ky為水平向和垂直向的滲透系數(shù)；Ψ為土水勢，參考文獻[3]將土體凍結和融化過程中的土水勢均定義為Ψ=Ψm+Y，其中Ψm為基質勢，Y為重力勢。

(3)

1.4 水熱聯(lián)系方程

上述公式(1)和公式(3)共有T，θu和θi三個未知量，只有兩個方程，因此需要補充一個聯(lián)系方程才能解方程組。聯(lián)系方程為未凍水含量與溫度變化的關系，如下所示：

θu=θ0|Ttrans|B|T|-B

(4)

式中：θ0為初始含水率Ttrans為冰水相變溫度；B為與土質相關的常數(shù)。

給定式(1)、式(3)和式(4)邊界條件和初始值就組成了路基工程水熱耦合的數(shù)學模型。

將未凍水含量θu對時間和空間求偏導得：

(5)

將水分遷移方程代入熱傳導方程并聯(lián)系上式可消去含冰率θi：

(6)

式(6)實現(xiàn)了方程組的解耦。將上式與水熱聯(lián)系方程(4)聯(lián)立即可求得溫度和水分的分布。

凍結區(qū)和未凍區(qū)之間有一個固態(tài)和液態(tài)分界界面，而且這個界面不固定。為了將該界面簡化為固定狀態(tài)，引入Heaviside階梯函數(shù)H(T)[15]利用COMSOL Multiphysics軟件內部自帶的表達式表示如下：

H(T)=flc2hs(T-Ttrans,dT)

(7)

式中：Ttrans為相變點溫度；dT為轉變溫度間隙取Ttrans=-0.3℃，dT=0.3℃，建立的函數(shù)變化曲線如圖1所示。

圖1 階梯函數(shù)曲線

根據(jù)式(7)，在整個求解區(qū)域內表示未凍水含量為：

θu=θ0|Ttrans|B|T|-B[1-H(T)]+θ0H(T)

(8)

容積熱容為：

C=Cf+(Cu-Cf)H(T)

(9)

式中：Cu為未凍土容積熱容量；Cf為已凍土容積熱容量。

滲透系數(shù)為：

k=kf+(ku-kf)H(T)

(10)

式中：ku為未凍土滲透系數(shù)；kf為已凍土滲透系數(shù)。

擴散系數(shù)為：

D=Df+(Du-Df)H(T)

(11)

式中：Du為未凍土水分擴散系數(shù)；Df為已凍土水分擴散系數(shù)。

參考以往研究成果[2,10]，定義阻抗系數(shù)I:

I=1010θi

(12)

那么滲透系數(shù)和水分擴散系數(shù)定義如下：

(13)

(14)

(15)

(16)

2 土體開放系統(tǒng)單向凍結試驗

為驗證土體開放系統(tǒng)中水熱耦合理論的有效性，選取路基填土制作圓柱形試樣，開展開放系統(tǒng)條件下單向凍結室內試驗。試樣直徑為38 mm，高度為76 mm。試樣初始含水率為20.2%，初始干密度為1.57 g/cm3。試樣制作完成后，用防水薄膜將試樣密封，僅保留底面以方便補水。試樣放置在圓柱形筒中，并安裝在如圖2所示的凍脹裝置中。該裝置能夠從上向下單向凍結試樣，最低溫度可達-40℃。此外，凍結過程中保持試樣底部恒溫在1℃左右，方便從外部給試樣補水。沿試樣高度方向每隔1.5 cm插入一個小型熱電阻，實時監(jiān)測試樣內部溫度變化。凍結試驗共進行8 h，試樣取出后，將土柱沿高度方向分層切片取樣烘干測含水率。試樣不同高度處的溫度隨凍結時間的變化及溫度分布如圖3所示?？梢钥吹皆嚇釉诤銣貎鼋Y的情況下，前2 h溫度劇烈變化，之后逐漸穩(wěn)定，凍結5 h后試樣全部達到0℃以下。并且隨著凍結時間的增加溫度分布呈現(xiàn)由非線性分布向線性分布過渡的趨勢，表明溫度梯度逐漸穩(wěn)定。

圖2 土體凍融循環(huán)試驗原理圖

3 水熱耦合理論的驗證

3.1 計算模型及參數(shù)

本次研究借助于第二節(jié)中開放系統(tǒng)土柱單向凍結試驗的溫度和水分變化結果驗證水熱耦合理論在開放系統(tǒng)季節(jié)凍土區(qū)應用的合理性。采用有限元軟件COMSOL Multiphysics根據(jù)圓柱土樣尺寸建立有限元模型如圖4所示，計算模型為二維軸對稱模型，其中高度為7.6 cm，半徑為1.9 cm。

圖3 凍結過程中不同時刻試樣的溫度分布

圖4 有限元計算模型

將試驗開始時室內溫度作為試樣的初始溫度，即T0=17.5℃。試樣初始含水率為20.2%，因此初始體積含水率設定為θ0=0.32。

試樣頂部即冷端(AD)溫度邊界條件按照試樣實際的溫度變化(冷卻溫度為-30℃時)進行設置。根據(jù)圖5可以發(fā)現(xiàn)試樣頂部(h=7.6 cm)溫度變化趨勢為先線性下降然后逐漸穩(wěn)定，為簡化計算過程，將試樣頂部溫度變化擬合為分段函數(shù)即式(17)，擬合結果如圖5所示。

(17)

圖5 試樣頂部溫度變化

AB為軸對稱邊界，CD為絕熱邊界，BC為暖端恒溫邊界，溫度為試驗過程中的實際溫度，Tb=2℃。對于水分邊界，AD和CD均為不透水邊界，AB為軸對稱邊界，BC為補水邊界，含水率恒定且為試樣的飽和體積含水率，即θs=0.42。參考徐學祖等[16]總結的亞黏土熱參數(shù)，綜合取值如表1所示[17]。那么可以得到土體的未凍水含量隨溫度的變化曲線如圖6所示。

表1 土體熱參數(shù)取值

圖6 未凍水含量隨溫度的變化曲線

COMSOL Multiphysics軟件自帶有系數(shù)型偏微分方程接口，可以將上述計算條件、參數(shù)及水熱耦合方程導入軟件中對方程重新編譯并進行計算。

3.2 結果與分析

根據(jù)室內試驗結果，試樣凍結持續(xù)5 h后，內部溫度全部達到零下并逐漸穩(wěn)定。因此以試樣凍結時間分別為1 h、2 h、4 h和5 h的溫度變化、凍結鋒面位置和水分分布對計算結果進行驗證。

(1)試樣溫度變化對比。試樣在凍結過程中的溫度變化云圖如圖7所示?？梢悦黠@看到隨著凍結時間的持續(xù)試樣凍結區(qū)逐步向下發(fā)展。數(shù)值計算時試樣底部保持恒溫2℃，因此圖7中凍結5 h后試樣的底部仍高于0℃。而試驗過程中雖然底板溫度也保持在2℃左右，但是凍融循環(huán)試驗系統(tǒng)中酒精循環(huán)液的冷卻效率明顯高于其保持恒溫的效率，因此出現(xiàn)5 h后試樣被完全凍結。

圖7 試樣凍結過程中的溫度變化云圖(單位：℃)

假定試樣凍結鋒面處土體的溫度為0℃，那么試樣在不同凍結時間計算凍結鋒面高度和實際凍結鋒面高度對比如圖8所示，其中試樣實際的凍結鋒面用白線標記。可以看到凍結持續(xù)1 h后凍結鋒面的計算高度和實際高度均為55 mm左右；凍結2 h后高度均為30 mm左右；凍結4 h后凍結鋒面的高度出現(xiàn)一定差異，實際高度為20 mm左右，而計算高度為15 mm左右；凍結5 h后室內試驗的試樣已經被完全凍結，而模擬結果的凍結鋒面的高度為5 mm左右。因此，在試樣完全凍結以前凍結鋒面的計算高度基本上和實際高度相差不大。

試樣不同高度處的計算溫度和實測溫度隨凍結時間的變化對比如圖9所示?？梢钥吹?，計算溫度和實測溫度變化趨勢一致，前2 h溫度迅速降低，凍結5 h后逐漸穩(wěn)定，且數(shù)值也比較接近。上述結果表明利用水熱耦合理論得到的開放系統(tǒng)土體凍結過程中的溫度計算結果是可信的。

圖8 試樣凍結過程中凍結鋒面的變化

圖9 試樣不同高度計算溫度和實測溫度變化對比

(2)試樣含水率變化對比。圖10為不同凍結時間試樣總體積含水率的計算值與實測值的對比。由于室內試驗中5 h后試樣即被完全凍結，而數(shù)值計算中溫度邊界條件的影響，土柱不會完全凍結。因此僅研究凍結時間分別為1 h、2 h和4 h時試樣內部的總體積含水率分布狀況。所述總體積含水率包含了未凍水含水率和冰晶融化為液態(tài)水后的體積含水率。從圖10中可以看到計算含水率與實測含水率僅在試樣頂部分布差異較大其余位置數(shù)值大小接近，整體的變化趨勢基本一致，均呈現(xiàn)試樣上部含水率高于下部并在中部發(fā)生突變，并且計算得到的含水率突變位置與實際的突變位置基本一致。將圖10與圖8對比，可以發(fā)現(xiàn)含水率突變位置即凍結鋒面位置，凍結鋒面以上為凍結區(qū)，凍結鋒面以下為未凍區(qū)。由于土體中的水分是由未凍區(qū)向凍結區(qū)遷移的，因此試樣上部含水率要高于下部。上述結果表明，利用水熱耦合理論計算土體的含水率變化是符合實際的，可以應用于季節(jié)凍土區(qū)受地下水遷移影響的巖土體的溫度和含水率計算中。

4 結 論

(1) 開放系統(tǒng)中土體的單向凍結試驗表明，在凍結過程中隨著凍結時間的增加，沿試樣高度方向的溫度分布呈現(xiàn)由非線性分布向線性分布過渡的趨勢。

(2) 室內試驗結果和水熱耦合模型的計算結果均顯示，土體在凍結期間未凍區(qū)水分在溫度梯度作用下向已凍區(qū)遷移，并在凍結鋒面處水分聚積，含水率發(fā)生突變。

圖10 不同凍結時間試樣的含水率分布

(3) 建立的開放系統(tǒng)條件下的土體水-熱耦合理論模型能夠較好地預測土體凍結過程中的溫濕度變化，可以應用于實際工程的溫度和含水率計算。