李 楊，王 清，王壇華

1.福建工程學院土木工程學院，福州 350014 2.吉林大學建設工程學院，長春 130026

?

凍土水熱耦合模型數(shù)值求解及結(jié)果檢驗

李 楊1,2，王 清2，王壇華2

1.福建工程學院土木工程學院，福州 350014 2.吉林大學建設工程學院，長春 130026

首先對作者所建立的基于多孔介質(zhì)理論的季節(jié)凍土水熱遷移耦合模型進行數(shù)值求解;對模型方程進行修正，并給出了模型方程中參數(shù)的確定方法。然后以長春——松原公路段土體為研究對象，對實際工程中凍結(jié)情況下水分遷移的情況進行預測;給定模型邊界條件對模型求解，將結(jié)果與野外實際監(jiān)測結(jié)果進行對比。溫度變化對比數(shù)據(jù)表明，模型可以較好地預測終值情況，而中間過程的誤差較大，但是趨勢基本一致。水分遷移方向及量的對比數(shù)據(jù)表明，模型計算結(jié)果要小于實測結(jié)果，但是整體上計算結(jié)果與實測結(jié)果的變化趨勢較一致，且同樣是和最終值吻合較好，誤差最小。結(jié)果表明，模型計算結(jié)果可較好地模擬參數(shù)最終值，但存在一定誤差。

多孔介質(zhì)；水熱遷移；耦合模型；溫度；含水率

0 引言

凍土是一種特殊的土體，它的固相由固體顆粒和冰共同組成。凍土可以視為多孔介質(zhì)，凍土中水分的遷移流動亦屬于多孔介質(zhì)流體流動，對凍土中水分遷移的研究也可以應用多孔介質(zhì)理論[1-4]。因此，凍土的水、熱遷移與成冰過程的本質(zhì)即為多孔多相介質(zhì)帶相變的固、液、氣、熱耦合問題。1957年，Philip和de Vries[5]基于多孔介質(zhì)中液態(tài)水黏性流動及熱平衡原理，提出了水熱耦合遷移模型，開創(chuàng)了土中水熱耦合研究的先河。此后Aboustit等[6]相繼對凍土中水熱遷移耦合模型進行了研究。國內(nèi)對凍土水熱耦合遷移問題的研究起步較晚，當前主要集中在研究正凍土尤其是飽和正凍土中水熱遷移問題等方面[7]，對非飽和季節(jié)凍土的水熱遷移耦合問題研究較少[8-11]，目前仍存在一些諸如不能很好地模擬自然界邊界條件、水分遷移的機制不確定及水分遷移過程的數(shù)值模擬與模型計算過于簡單等問題。本文對筆者建立的一個基于多孔介質(zhì)理論的水熱遷移耦合模型[12-13]進行求解，并將結(jié)果與實例進行對比，從而驗證模型的可靠性，以使模型能更好地模擬實際情況，解決實際問題。

1 模型求解基礎

1.1 原始模型方程

所求解的原始模型方程[13]如下：

式中：t為時間；T為溫度；ρd為土粒干密度；e為孔隙比；ρs為土粒密度；ρl為水的密度；θl,θi,θ分別為液態(tài)含水率、含冰率以及總含水率；C為骨架熱容；hi,l為水的凍結(jié)(或冰融解)潛熱；λl,λi,λs分別為液態(tài)水、冰及土骨架導熱系數(shù)；cl,ci分別為液態(tài)水和冰的比熱；Dl為液態(tài)水擴散系數(shù)；Kl為液態(tài)水滲透系數(shù)。

1.2 修正模型方程

為了方便對模型進行求解，對模型進行如下修正。根據(jù)控制方程式(1)和(2)，考慮到由于溫度變化造成水分的相變和遷移，引入階躍函數(shù)H(x)(Heaviside函數(shù))作為相變指示函數(shù)f(T)，有

則總含水率為

θ=θlf(T-273)+θif(273-T) ；

水分控制方程(1)可改寫為

Dleρs2θl+ρl

溫度傳遞方程(2)可改寫為

λlf(T-273)2T+λif(273-T)2T+

將式(6)兩邊同時乘以hi,l再代入(7)得

λlf(T-273)2T+λif(273-T)+

hi,lρlKl(θl)+λs2T+

根據(jù)未凍含水率和溫度的關(guān)系式θl=θl(t,T)，得

式(9)代入式(8)得

[λlf(T-273)+λif(273-T)+λs+

f(T-273)+hi,lρlk

代入階躍函數(shù)H(x)，由式(7)得

(273-T)-clρl

代入階躍函數(shù)H(x)，由式(8)得

hi,lρlKf(T-273)=λlf(T-273)2T+

λif(273-T)2T+λs2T+

clρl

式(5)、(10)、(11)、(12)即為修正的模型控制方程，對方程采用有限元法進行二維求解，即可得到數(shù)值解。

1.3 模型參數(shù)確定

為了保持單位的一致性，文中所有參數(shù)一律采用國際單位制。部分參數(shù)采用在長春——松原一級公路線段采取的土樣進行相關(guān)試驗及觀測后得到的數(shù)值。

1)已知參數(shù)

λl=0.574 5 W/(m·K)，λi=2.1 W/(m·K)，cl(288 K)=1.002 J/(kg·K)，ci=0.505 J/(kg·K)，hi,l=334 kJ/kg。

2)實驗測定參數(shù)

液態(tài)水擴散系數(shù)與體積含水率的關(guān)系式為

相關(guān)度為0.944 8。

土體導熱系數(shù)與溫度的關(guān)系擬合式為

λs=-0.006 6T+3.198 9。

未凍含水率與溫度的關(guān)系擬合式為

θl= 4.1e-64T25.696。

導水系數(shù)與未凍含水率的擬合關(guān)系式為

K=Cw(1.205 4e+4θl10.5)，

其中，Cw=1 500 kJ/(m3·K)。

當溫度為291 K時，骨架比熱cs=830 J/(kg·K)，當溫度為263 K時，cs=763 J/(kg·K)。

C可視為體積含水率和體積含冰率的函數(shù)：

2 模型求解

2.1 邊界條件

模型的求解以長春地區(qū)的參數(shù)為基礎。長春地區(qū)冬季嚴寒、漫長，凍結(jié)期長達5個多月。標準凍深約1.7 m，為典型的溫帶氣候區(qū)。選取長春——松原公路段土樣的相關(guān)參數(shù)，利用所建立的模型對其進行求解，模擬長春地區(qū)季節(jié)凍土的凍結(jié)發(fā)展情況，研究土體中不同深度、不同溫度下的水分遷移過程。

參考長松線公路野外監(jiān)測結(jié)果，將模型的深度下限定為3.0 m，寬度取1.0 m。因為所監(jiān)測地點為公路地段，地面有50 cm的路面層，監(jiān)測數(shù)據(jù)為從地表以下0.5 m處開始，所以計算模型上部邊界選取深度為0.5 m處。由于本區(qū)地下水埋深在7.0 m左右，埋深較大，毛細水無法為凍結(jié)區(qū)補給水分，且熱量和孔隙水均沿垂直向運移，水平向的熱和孔隙水交換量很小，所以可以將水熱遷移視為一維問題。模型初始溫度模型示意圖如圖1所示。模型條件如下。

圖1 土柱模型示意圖Fig.1 Chart of the soil column model

1)初始條件

溫度T(t=0,y)=10.3 ℃=283.3 K；含冰率θi=0；總含水率θ=2.571y4-11.285y3+10.584y2+3.810y+19.422。

2) 邊界條件

模型頂端溫度屬于第一類邊界條件，但由于溫度隨時間變化，根據(jù)現(xiàn)場實際監(jiān)測結(jié)果進行擬合得到

T(t,y=0.5)=264.3+3e-6t。

底端溫度屬于第二類邊界條件：

總含水率θ、含冰率θi在模型兩端均為第二類邊界條件：

2.2 模型的解

模型采用有限元進行求解，模擬的過程包括從剛?cè)攵瑫r的不穩(wěn)定凍結(jié)階段一直到氣溫回升之后的雙向融化階段(即2004——2005年)，共5個月，所得結(jié)果如下。

2.2.1 凍結(jié)過程中的溫度變化

圖2為凍結(jié)1～5月的土柱溫度分布曲線，圖3為凍結(jié)第5月的土柱溫度分布和傳遞方向。由圖2可見：1)凍結(jié)第1月，土體溫度分布比較均勻，上部溫度較低，下部溫度較高，此時凍結(jié)鋒面位于頂部。這是因為現(xiàn)在的趨勢為在土柱中部分開，下部土體溫度向上部傳遞，上部土體溫度向下部傳遞。2)凍結(jié)第2月，土體上部負溫區(qū)(T<237 K)開始增大，說明溫度向下部傳遞。而此時凍結(jié)鋒面已經(jīng)明顯下移。3)到第3月時，土體負溫部分繼續(xù)下移，在土體最上部以及中間部位開始出現(xiàn)升溫的現(xiàn)象，但并不明顯。4)到達第4月時，土柱負溫區(qū)開始縮小，并且在土柱上部開始有正溫區(qū)的出現(xiàn)，此時凍結(jié)鋒面以上的溫度傳遞方向為負溫向上傳遞。凍結(jié)鋒面以下土體溫度繼續(xù)向下傳遞。5)到第5月時，地表溫度早已達到正溫，從土柱中也可以看出，正溫區(qū)明顯增多，負溫區(qū)變小，此時地表出現(xiàn)溫度不穩(wěn)定區(qū)。從圖3亦可以看出,此時土柱中溫度傳遞方向整體向下,但靠近地表處溫度傳遞方向略顯凌亂,與圖2顯示了相同的趨勢。

2.2.2 凍結(jié)過程中的含水率變化

圖4為凍結(jié)1～5月的土柱含水率分布曲線，圖5為為凍結(jié)第5月的土柱含水率分布和傳遞方向。由圖4可見：1)在第1月，土柱中以1.0 m深度左右為分界線，上部土體中的水分向下遷移，下部土體中的水分向上遷移。2)到第3月時，可以看出水分遷移仍具有第1月時的趨勢，且此時的1.0 m上部的含水率數(shù)值比第1月時有所增大。1.0 m以下部分含水率均有所減小。3)到第5月時此趨勢更加明顯，地表含水率的數(shù)值比凍結(jié)第1月時明顯增大，1.0 m下部的也繼續(xù)減少。由圖5也可以看出,土柱中1.0 m以上水分向上遷移,而1.0 m以下有一個向下遷移的趨勢。這是因為隨著溫度的降低，土體中負溫開始從地表向下傳遞，同時土體由地表向下凍結(jié)，當上部土體凍結(jié)成冰后，其下部土體中的水分必然向上遷移對其補充。本區(qū)地下水埋深在7.0 m以下，所以毛細水無法供給水分遷移，因此導致1.0 m下部的土體含水率減小。

3 解的檢驗

3.1 計算溫度與監(jiān)測值對比

為了對求解結(jié)果進行檢驗，將計算值與監(jiān)測值進行對比。文中所采取的地溫監(jiān)測數(shù)據(jù)為與計算模型采用的土樣一致的長松線2004——2005年度的地溫監(jiān)測結(jié)果。地溫監(jiān)測采取鉆孔分段埋設探頭，根據(jù)需要按規(guī)定的時間進行讀數(shù)，精度0.5 ℃。

凍結(jié)第1、3、5月時計算的土體溫度與實測的土體溫度對比見圖6。從圖6可以看出，整體上模型計算值與實測值的趨勢是一致的，但是存在誤差，且隨著凍結(jié)時間增長，誤差值有減小的趨勢，模型計算值更加接近實測值。

凍結(jié)第1、3月時皆是計算值中土層深度存在一個分界區(qū)，此區(qū)域上部模型計算值比實測值大，下部計算值比實測值小。在凍結(jié)第1月時，總體上以深度為2.0 m處為分界線，2.0 m以上計算結(jié)果比實測值大，2.0 m以下計算結(jié)果比實測值小，最大誤差值達2 ℃左右。第3月時分界區(qū)域有所上移，但趨勢基本一致,且誤差有比凍結(jié)第1月時減小的趨勢。

土體凍結(jié)第5月，即計算與觀測的最終結(jié)果時，從圖6中可以看出，模型計算結(jié)果與實測結(jié)果吻合得較好，變化趨勢基本一致，誤差只有1.0 ℃左右。從對地溫要求的角度來說，這個誤差范圍是可以接受的，因為更重要的是了解溫度變化的趨勢，因此誤差滿足要求。

圖2 凍結(jié)1～5月的土柱溫度分布曲線Fig.2 Temperature curve of the soil column frozen for 1-5 months

圖3 凍結(jié)第5月的土柱溫度分布和傳遞方向Fig.3 Temperature gradient distribution and transmission direction of the soil column frozen for five months

圖4 凍結(jié)1～5月的土柱含水率分布曲線Fig.4 Moisture content curve of the soil column frozen for 1-5 months

圖5 凍結(jié)第5月的土柱含水率分布和傳遞方向Fig.5 Moisture content gradient distribution and transmission direction of the soil column frozen for five months

從以上3個不同時間段計算結(jié)果與實測結(jié)果的對比結(jié)果可以得出，所建立的模型可以較好地預測終值情況，雖然中間的情況誤差較大，但是趨勢基本一致。

可能導致出現(xiàn)這種現(xiàn)象的原因很多。例如：監(jiān)測深度的密度不夠大，以致于有些可能存在的溫度突變點沒有監(jiān)測出來；模型用于計算的實測溫度變化曲線為對現(xiàn)有的幾個溫度監(jiān)測數(shù)據(jù)進行擬合的結(jié)果，必然存在一定的誤差；土體的各種參數(shù)都是取樣后室內(nèi)測定的，存在一定誤差；在現(xiàn)實中，天氣等各方面因素的變化都可能導致氣溫的不穩(wěn)定性，氣溫的變化必然導致土體溫度的變化；而且路基上道路行車等給路基帶來的荷載效應等都會導致計算結(jié)果與實測結(jié)果有差異。因此，為使模型能更加符合千變?nèi)f化的實際情況，今后應該對這些方面做進一步的研究。

3.2 計算土體含水率與實測值對比

文中所采取的實測值為與計算模型采用的土樣相同的長松線2004——2005年度的含水率監(jiān)測結(jié)果。

土體的計算含水率值與實測值在凍結(jié)第1、3、5月時的比較情況見圖7。由圖7可知，3個凍結(jié)時間段基本上都是計算結(jié)果小于實測結(jié)果，但是整體上計算結(jié)果與實測結(jié)果的變化趨勢還是比較一致的；和溫度計算結(jié)果一樣，同樣是土層中部深度處模擬較好，上部和下部誤差較大，而且同樣模型的計算結(jié)果和最終值吻合得較好，誤差最小，平均只有2%左右，這說明模型計算結(jié)果的終值能較好地模擬實際情況。

圖6 凍結(jié)第1、3、5月土柱溫度實測與計算對比圖Fig.6 Temperature comparison chart of measurement and calculation frozen for one, three and five months

圖7 凍結(jié)第1、3、5月含水率實測與計算對比圖Fig.7 Moisture content comparison chart of measurement and calculation frozen for one, three and five months

4 結(jié)論與討論

對作者建立的水熱遷移耦合模型進行修正，以長春——松原公路路段土樣為基礎，確定了模型參數(shù)，在一定假設的基礎上，進行了數(shù)值求解，得到了土體凍結(jié)過程中溫度及含水率在不同深度以及不同時間上的分布、變化情況。并將求解結(jié)果與長春——松原公路路段實際監(jiān)測結(jié)果進行了對比。

根據(jù)對比結(jié)果，對于溫度和含水率，模型的計算結(jié)果都與實測結(jié)果的趨勢基本一致，但存在一定誤差，整體上位于精度范圍內(nèi)。誤差出現(xiàn)的原因可能如下：1)模型的建立是為了對現(xiàn)實情況進行模擬，但自然界是千變?nèi)f化的，不可能對其所有的參數(shù)進行量化，且模型的建立需要考慮的因素甚多，在模型建立的過程中，為了方便計算、求解，必然會對實際情況進行一些簡單化處理，也就是預先設定一系列的假設。本次所求解的模型在建立時就進行了一系列的假設，這必然會給模型的計算結(jié)果帶來誤差。2)模型的求解是建立在實測的土體參數(shù)基礎上的，因此參數(shù)測定的準確與否對模型解的正確性也起著決定性的作用。凍土作為一種復雜的多相體，涉及到水熱的對流、傳導、相的轉(zhuǎn)變等一系列物理化學過程，而能代表這些變化的熱交換參數(shù)和水力參數(shù)測定的準確性卻至今是個難點，而這些測定結(jié)果又對最終的計算結(jié)果影響較大，因此也會帶來較大的誤差。且土體隨著深度的變化，各種參數(shù)也會有所變化，本次研究中為了方便模型的求解，對土體各個深度都采用了相同的參數(shù)，這樣的處理必然也會給計算結(jié)果帶來誤差。3)模型的計算結(jié)果是與實際監(jiān)測值進行對比，而實際監(jiān)測值受監(jiān)測方法、儀器精度等各種因素的影響，本身也存在一定誤差。4)現(xiàn)實的自然界是復雜多變的，道路上覆積雪的融化向下滲流會對含水率的分布造成很大影響，而且作為路基，土體中的水分遷移過程將受車輛等外部動荷載作用，同時土體變形(壓密、凍脹)造成的應力場對水分遷移也有影響。而這些在模型中都是無法反映的，因此計算結(jié)果與實測結(jié)果會有一定的誤差。

多種因素可導致模型計算結(jié)果與實測結(jié)果有出入，如何在模型中將各種影響因素考慮進去以縮小誤差，將結(jié)果控制在合理的區(qū)間內(nèi)需要進一步研究。例如：可以在模型中考慮積雪融化入滲對水分遷移的影響；利用更合理、準確的方法對土體的各種參數(shù)進行測定；可以將由凍脹產(chǎn)生的力引入到模型中建立水熱力耦合模型等。這些都是將來可以繼續(xù)研究的問題。

[1] 朱志武，寧建國，馬巍. 土體凍融過程中水、熱、力三場耦合本構(gòu)問題及數(shù)值分析[J]. 工程力學，2007，24(5)：138-144. Zhu Zhiwu, Ning Jianguo, Ma Wei. Constitutive Model and Numerical Analysis for the Coupled Problem of Water, Temperature and Stress Fields in the Process of Soil Freeze-Thaw[J]. Engineering Mechanics，2007, 24(5)：138-144.

[2] 羅煥炎，陳雨孫. 地下水運動的數(shù)值模擬[M]. 北京：中國建筑工業(yè)出版社，1988. Luo Huanyan, Chen Yusun. Numerical Simulation of Groundwater Movement[M]. Beijing: China Building Industry Press,1998.

[3] Bear J, Bechmat Y. Introduction to Modeling of Transport Phenomena in Porous Media[M]. [S.l.]:Kluwer Academic Publishers, 1991.

[4] 林瑞泰. 多孔介質(zhì)傳質(zhì)傳熱引論[M]. 北京：科學出版社，1995. Lin Ruitai. Heat and Mass Transfer in Porous Media Introduction[M]. Beijing: Science Press,1995.

[5] Philip J R, De Veries D A. Moisture Movement in Porous Material Under Temperature Gradient[J]. Trans Am Geophys Union, 1957, 38:222-232.

[6] Aboustit B L, Advani S H, Lee J K, et al. Finite Element Evaluation of Thermo-Elastic Consolidation[J]. Proc US Symp Rock Mech，1978, 23:587-595.

[7] 周家作, 李東慶, 房建宏, 等. 開放系統(tǒng)下飽和正凍土熱質(zhì)遷移的數(shù)值分析[J]. 冰川凍土, 2011, 33(4):791-795. Zhou Jiazuo, Li Dongqing, Fang Jianhong, et al. Numerical Analysis of Heat and Mass Transfers in Saturated Freezing Soil in an Open System[J]. Journal of Glaciology and Geocryology, 2011, 33(4):791-795.

[8] 陳飛熊, 李寧, 程國棟. 飽和正凍土多孔多相介質(zhì)的理論構(gòu)架[J]. 巖土工程學報, 2002, 24(2): 213-217. Chen Feixiong, Li Ning, Cheng Guodong. Theoretical Framework Freezing Soil Saturated Porous Media Multiphase[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(2): 213-217.

[9] 白冰, 劉大鵬. 非飽和介質(zhì)中熱能傳輸及水分遷移的數(shù)值積分解[J]. 巖土力學, 2006, 27(12): 2085-2089. Bai Bing, Liu Dapeng. Numerical Integral Solutions of Heat Transfer and Moisture Transport in Unsaturated Porous Media[J]. Rock and Soil Mechanics, 2006, 27(12): 2085-2089.

[10] 吳曉玲, 向小華, 王船海, 等. 季節(jié)凍土區(qū)融雪凍土水熱耦合模型研究[J]. 水文，2012(5):12-15. Wu Xiaoling, Xiang Xiaohua, Wang Chuanhai, et al. Numerical Modeling for Coupled Snowmelt and Frozen Soil in Seasonally Frozen Soil Region[J]. Journal of China Hydrology, 2012(5):12-15.

[11] 宿曉萍，王清，王文華，等. 季節(jié)凍土區(qū)鹽漬土環(huán)境下混凝土抗凍耐久性機理[J]. 吉林大學學報:地球科學版，2014，44(4)：1244-1253. Su Xiaoping, Wang Qing, Wang Wenhua, et al. Frost Resistance and Durability Mechanism of Concrete Under Saline-Alkali Condition in Seasonal Frozen Soil Area[J]. Journal of Jilin University：Earth Science Edition, 2014, 44(4): 1244-1253.

[12] 李楊.基于多孔介質(zhì)理論的季節(jié)凍土水熱遷移耦合模型推導[J]. 河北工程大學學報:自然科學版，2012(3): 11-14. Li Yang. Derivation of Water and Thermal Migration Coupled Model Based on Theory of Porous Media in Seasonal Frozen Soil[J]. Journal of Hebei University of Engineering:Natural Science Edition, 2012(3): 11-14.

[13] 李楊. 季節(jié)凍土水分遷移模型研究[D]. 長春: 吉林大學, 2008. Li Yang. A Study on the Moisture Content Migration Model of Seasonal Frozen Soil[D]. Changchun: Jilin University, 2008.

Numerical Solution and Test of Results for a Hydrothermal Coupled Model About Frozen Soil

Li Yang1,2, Wang Qing2, Wang Tanhua2

1.DepartmentofCivilEngineering,FujianUniversityofTechnology,Fuzhou350014,China
2.CollegeofConstructionEngineering,JilinUniversity,Changchun130026,China

A hydrothermal coupled model established previously by the author based on the theory of porous media seasonally frozen soil was solved numerically in this article. First, the model equations were corrected, and the method to determine parameters in equations was given. Then Changchun-Songyuan highway soil was chosen as the object for the study, and the case of frozen moisture migration can be predicted in the actual project. Model was solved under the given circumstances the model boundary conditions, and the model solution results were compared with actual field monitoring results. The comparative data of temperature gradient show that the model can predict the final value of the parameter, while the middle of the process error is large, but the same trend. The comparative data of direction and the amount of moisture migration indicates that the model results to be less than the measured results, but overall the same trend. And the same with temperature gradient that the final values agree well with the smallest error, an average of about 2%. The results show that the model can simulate calculated parameters for the final values, but there is an error, and the error is within an acceptable range.

porous media；hydrothermal transfer；coupling model；temperature；moisture content

10.13278/j.cnki.jjuese.201501202.

2014-04-27

國家自然科學基金項目(40672180,41372267)；高等學校博士學科點專項科研基金項目(20120061110054)；福建工程學院科研啟動基金項目(GY-Z12077)

李楊(1981——)，女，工程師，博士，主要從事工程地質(zhì)、巖土工程，地質(zhì)災害防治等方面的研究，E-mail:liyang513@foxmail.com。

10.13278/j.cnki.jjuese.201501202

TU443

A

李楊，王清，王壇華. 凍土水熱耦合模型數(shù)值求解及結(jié)果檢驗.吉林大學學報:地球科學版,2015,45(1):207-213.

Li Yang, Wang Qing, Wang Tanhua. Numerical Solution and Test of Results for a Hydrothermal Coupled Model About Frozen Soil.Journal of Jilin University:Earth Science Edition,2015,45(1):207-213.doi:10.13278/j.cnki.jjuese.201501202.