寒區土體一維水-熱-力耦合模型與數值分析

魏道凱，荊 皓，陳 琦，寇海磊

（1.山東交通職業學院公路與建筑系，山東濰坊 261206；2.中國海洋大學工程學院，山東青島 266100）

引言

凍土在我國國土面積中占據著極其廣闊的一部分，其中多年凍土約占1/4。寒區凍土的凍脹融沉特性使結構物產生不均勻沉降、結構斷裂等現象，工程災害頻頻發生，對上覆結構物的穩定帶來了極大的挑戰。凍土凍融作用的機理是當前凍土研究中最受關注的問題之一［1-5］。

在外界溫度環境作用下，土體內液態水會因為土體內部溫度低于土體的凍結溫度引起物相改變，水變成冰，冰的體積較水的體積大，顆粒間的孔隙減小，進而實現土的凍脹。同時，冰的膠結作用使土顆粒間應力增大。故凍土凍融為溫度場、水分場和應力場三場耦合的結果。對于三場耦合的研究一直是國內外研究學者關注的重點。Harlan［6］首次將達西定律應用于凍土水分場分析中，建立水熱耦合方程。Mu等［7］首次提出凍土中水熱兩場的變化會引起土體內部應力場的變化。近些年關于凍土凍融機理的研究都是在此基礎上進行的。Zhou等［8］基于人工凍結試驗將土體凍脹引入水熱力耦合模型中，探究了水熱作用下土體變形特性?；谒疅崃θ龍鲴詈夏Ｐ?，Qi等［9］對凍融環境下膨脹土邊坡的穩定性進行了分析。周志［10］通過水熱耦合模型對凍融環境下粉質黏土內溫度場以及水分場進行了研究。冉洪伍等［11］對水動力學模型、剛性冰模型以及熱力學模型3種耦合模型的應用以及準確性展開了研究。但是已有研究多針對于飽和土體，對于非飽和凍土在水熱兩場作用下應力場以及凍脹變形的研究較少。

文中基于質量守恒定律以及能量守恒定律推導水分場控制方程及溫度場控制方程，由水分場和溫度場作用產生土體內部應力變化，實現土體的凍脹，建立與應力場的耦合。以此建立一維土柱數值模型，探究三場耦合的可行性。對凍融條件下土柱溫度變化、水分遷移以及應力分布進行了研究。

1 水-熱-力三場耦合數學模型

1.1 基本假設

在保證模型準確性的基礎上忽略部分情況，充分考慮凍土中水熱力因素，簡化模型建立過程，做出以下基本假設：

（1）在土體中，水分遷移主要以液態形式進行遷移，在模型中忽略氣態水的遷移；

（2）土體中水分的遷移服從廣義達西定律；

（3）土體各向同性，均勻連續；

（4）計算過程無溫度損失，凍土內含水量以及溫度處于平衡狀態。

1.2 水分場控制方程

根據質量守恒定律建立水分場控制方程。二維土體微單元水分遷移模型如圖1所示。

圖1 土體微單元水分遷移示意圖Fig.1 Schematic diagram of water migration in soil micro-units

式中：Δmx為單位時間土單元中x方向水分質量變化；vx為x方向水流通量；Δmz為單位時間土單元中z方向水分質量變化；vz為z方向水流通量；ρl為水的密度。

則單位時間內土體微單元的水分增量為：

單位時間內土體微單元水分質量變化有可從未凍水體積含量θl以及冰體積含量θi的變化進行分析。

聯立方程可得：

式中ρi為冰的密度。

水分遷移服從廣義達西定律，則水流通量v：

式中：k（θl）為土體的滲透系數；ψm代表基質勢；z代表重力勢。

最終得到水分場控制方程：

1.3 溫度場控制方程

基于能量守恒定律建立土柱內溫度場控制方程。

圖2表示土體微單元內的熱量遷移，與水分場類似，在單位時間內，土體微單元的能量變化ΔQ為：

熱通量q又為導熱系數λ與溫度T的函數：

當溫度降至土體凍結溫度時，土顆粒間液態水會發生物相改變以釋放能量［12-13］，因此土體微單元內能量變化又是由土體熱容C以及相變潛熱L決定的。

聯立方程建立溫度場控制方程：

水分場、溫度場控制方程中共有θl、θi以及T這3個未知函數因此需要還需一個聯系方程實現水熱兩場的耦合。白青波等［14］提出固液比B（T）的概念，將其作為聯系方程：

式中：Tf為土體的凍結溫度；b為固液比系數，砂土通常取值0.61，粉土取值0.47，黏土取值0.56。

文中水分場與溫度場之間屬強耦合，彼此相互影響。方程（15）以土體凍結溫度Tf為分界點可分為2個階段：當溫度超過土體的凍結溫度時，土體內不會發生冰水相變過程，此時溫度變化與土體的熱容以及導熱系數有關；當溫度低于土體的凍結溫度時，土體內水開始相變成冰。一方面，相應位置含水量減少，基質勢發揮作用，土體內水分開始遷移，由液態水含量高的位置遷移至液態水含量低的位置。另一方面，冰水相變所產生的潛熱也開始影響溫度的變化速率，水分場與溫度場之間的耦合開始發揮作用。

1.4 應力場控制方程

土體在凍結過程中，水分場重分布以及水相變成冰產生的凍脹造成土體內部應力重分布［15］。由水熱兩場耦合產生的體積應變作為土體內部的應力場，從而建立應力場控制方程。

圖2 土體微單元熱通量示意圖Fig.2 Schematic diagram of heat flux in soil micro-unit

靜力平衡方程：

幾何方程：

物理方程：

在凍結過程中，凍結區內的水分產生了相變，水相變成冰體積相應增大9%。與此同時，在基質勢的作用下，未凍結區的液態水向凍結區移動，則土體的應變可以用遷移的液態水含量來表示，即溫度作用下由于水分遷移產生的應變εvf為：

式中：θ0為初始含水量；Δθ為水分遷移量；θu為未凍水含量；n為初始孔隙比。

需要說明的是，土體凍融是一個非常復雜的物理過程，溫度場以及水分場變化影響應力的重分布。反之，應力重分布又影響著土體溫度以及體積含水量的變化［16］。在溫度場作用下，凍結區的水分發生相變導致凍結區液態水含量降低，進而在基質勢作用下未凍結區液態水遷移至凍結區，其填充土體內孔隙并產生進一步相變，水相變成冰體積增大，土體內部結構產生應力。在應力作用下土體發生應變，進而造成土體的凍脹。

2 基于Comsol數值模型構建

本次數值模型參數取自涂志斌［17］一維土柱室內模型試驗。參考試驗用土為膨脹土，土體參數如表1所示。試驗采用如圖3所示裝置，試驗開始前將土體分層填充至裝置中，最終形成直徑14 cm、高40 cm的土柱，采用直線位移傳感器測量試驗過程中土體頂部的凍脹量，采用土壤溫濕度一體傳感器測量距離土柱頂部5 cm位置處的溫度及含水量變化。對土柱頂部進行密封處理，將土柱置于多功能環境箱內，設置環境溫度為10℃并維持24 h，以確保土樣中的溫度均勻分布，然后在試驗裝置外側及底部包裹保溫棉以達到單向凍結的目的。最后，將試驗裝置放于多功能環境箱內，控制環境溫度為-10℃并維持144 h。

表1 試驗用土基本物理參數表［17］Table 1 Table of basic physical parameters of test soil

由于在土柱底板以及側邊施加保溫措施，本次數值模擬假定土柱四周及底部為絕熱狀態且完全側限，土柱與外界之間的換熱僅由頂端發生。由于溫度，水分只沿軸向傳遞，且只發生軸向變形，因此，本次數值模擬可簡化為一維土柱水-熱-力耦合模型。

數值模擬基于有限元數值軟件Comsol。為簡便計算，定義凍土相對飽和度S為：

式中：θu為未凍水含率；θs為飽和含水率，選擇S為變量代替θ進行求解。

將上述水分場、溫度場控制方程編寫成圖4（a）所示格式進行求解。其中，水分場控制方程為：

溫度場控制方程為：

對于應力場采用Comsol內置固體力學模塊進行求解。水分場邊界條件設為零通量，溫度場邊界條件選用狄式邊界條件。狄式邊界條件可以是隨時間變化的變量，也可以是恒定值，方程形式如圖4（b）所示，初始溫度為10℃，上邊界溫度為-10℃，下邊界設置為絕熱。模型計算中所需參數如表2所示。其中，土柱密度、初始體積含水量、彈性模量均根據參考試驗取值，土體的導熱系數以及熱容參數根據《凍土地區建筑地基基礎設計規范》（JGJ 118-2011）［18］建議值選取。為簡化計算，對土柱進行對稱建模，沿土樣高度劃分為450個域單元以及118個邊界單元，有限元數值模型如圖5所示。

圖3 監測點布置示意圖［17］（單位：cm）Fig.3 Layout diagram of monitoring points（Unit：cm）

圖4 Comsol內置假設方程形式Fig.4 Equation form in Comsol

圖5 有限元數值模型（單位：m）Fig.5 Finite element numerical model（Unit：m）

表2 數值模型參數表Table 2 Table of numerical model parameters

3 結果與討論

3.1 溫度場

圖6為距離土柱頂部5 cm位置處溫度隨凍結時間的變化。分析可知，溫度變化曲線可分為2個階段：在凍結時間12 h內，該位置處溫度下降速度較快，由10℃迅速下降至土體的凍結溫度（-1.33℃），這是由于土柱頂部位置溫度梯度較大，在此階段無冰晶出現，未發生冰水相變。待溫度達到凍結溫度后，冰晶成核迅速生長，土體內自由水發生相變并迅速凍結；在凍結時間超過12 h后，由于水相變吸熱，溫度下降速度變緩，在遞降階段中，土體中的結合水也開始發生凍結。數據顯示，溫度場數值模擬結果與實測結果具有較好的吻合度，驗證了溫度場控制方程以及耦合的準確性。

圖7顯示了在冷端溫度為-10℃時，不同制冷時間下一維土柱的溫度場有限元數值模擬云圖。隨著制冷時間的增加，凍結面逐漸向下傳遞，在凍結12 h時，凍結面距土柱底部0.35 m處，與圖6實測結果有較好的吻合度，在凍結144 h時，凍結面距土柱底部0.24 m。圖8表示不同凍結時間下溫度沿土柱高度分布。由圖可知，凍結區與未凍結區之間形成了凍結鋒面，且下移速度在開始的1～48 h內較快，這是由于在凍結初始階段土

圖6 距離土柱頂部5 cm處溫度隨凍結時間變化Fig.6 Temperature variation with freezing time at 5 cm from the top of soil column

圖7 溫度場云圖Fig.7 The cloud map of temperature field

內溫度變化幅度較大。由圖可以看出土柱內的溫度分布曲線中存在明顯的拐點，這是在溫度降至凍土的凍結溫度時，凍土內液態水在該溫度處發生物相改變，凍結區與未凍結區的溫度分布出現了明顯的差異。在凍結后期，凍結深度變化速度變緩。在凍結144 h后，溫度沿土柱高度近似呈線性分布。

3.2 水分場

圖9為距離土柱頂部5 cm位置處液態水含量隨凍結時間的變化。由圖可知，該位置處液態水含量隨凍結時間不斷降低。在凍結時間12 h內，液態水含量幾乎無變化，這是由于該位置處的溫度高于土體的凍結溫度，液態水未發生物相改變。當凍結時間超過12 h后，溫度降低至土體的凍結溫度，液態水發生物相改變，液態水含量降低，且前期降低速率較快。與此同時，在基質勢作用下液態水開始遷移，由土柱下部遷移至土柱上部。該位置處液態水在相變成冰的同時，又受到來自土柱下部液態水的補充，故后期液態水含量降低速率逐漸變緩。水分場數值模擬結果與實測結果具有較好的吻合度，并證明了水分場控制方程以及水-熱耦合的準確性。

圖10表示在不同凍結時間下土柱內含水量分布的數值模擬結果。由圖分析，在土柱凍結過程中，凍結區含水量明顯增加，在凍結區與未凍結區處含水量差異較大，表明未凍結區處的液態水不斷向凍結鋒面遷移。這是由于凍結區液態水在低溫環境下發生相變成冰，導致液態水含量減少，在水分場控制方程中，液態水會在基質勢的作用下由含量高的地方遷移至含量低的地方，這就造成在凍結鋒面處形成鮮明的S形曲線［19-20］。凍結鋒面在一維凍結的條件下隨凍結時間逐漸下移，這與溫度場的變化規律保持一致。

圖8 不同凍結時間土柱溫度分布圖Fig.8 Temperature distribution map of soil column with different freezing time

圖9 距離土柱頂部5 cm處液態水含量隨凍結時間變化Fig.9 Water content variation with freezing time at 5 cm from the top of soil column

圖10 不同凍結時間含水量分布圖Fig.10 Distribution map of water content at different freezing time

3.3 凍脹位移

圖11表示土體凍脹量隨凍結時間的變化關系。為便于比較，將考慮多場耦合以及不考慮多場耦合凍脹量計算結果同樣列于圖11中。分析可知，在凍結30 h前，土體凍脹變形發展較快。這是因為冷端溫度傳遞速度快，上部土體孔隙被冰迅速填充，進而發生凍脹。在凍結后期，凍脹變形速度變緩，最終凍脹量為5 mm。進一步分析可得，考慮多場耦合計算出數值模擬結果（即以式（19）計算的結果）與模型試驗結果所表現的凍脹量變化趨勢相似，且最終凍脹量均為5 mm，如圖12所示。然而，未考慮多場耦合計算出的數值模擬結果雖然變化趨勢與實測值相似，但計算出的凍脹量明顯小于實測值。當不考慮多場耦合時，即溫度場以及水分場對應力場無影響，此時的凍脹應變εvf為：

此時，凍脹應變僅僅為初始含水量條件下土體內水相變成冰所產生的應變，不包含溫度場以及水分場作用下的水分遷移所導致的應變增加，故未考慮多場耦合所計算出的數值模擬結果遠遠小于實測值。該對比結果進一步驗證了文中水熱力三場耦合的準確性以及可靠性。

圖11 凍脹量隨凍結時間變化Fig.11 Relationship between the frost heave changes and freezing time

4 結論

圖12 凍結144 h后土柱凍脹位移云圖Fig.12 Frost heave displacement cloud map of soil column after freezing for 144 h

文中基于質量守恒定律以及能量守恒定律對凍土水分場以及溫度場控制方程進行系統的推導，并結合水熱條件下凍土內部應力變化建立水熱力三場耦合數學模型，通過建立一維非飽和土柱數值模型可得結論：

（1）在凍結過程中，凍結面不斷下移，初始凍結時間下移速度較快，后期趨于平穩，最終土柱內溫度隨高度線性分布；

（2）土柱內水分在凍結過程中發生遷移，在凍結面處形成明顯的S型曲線，未凍結區水分不斷向凍結區遷移；

（3）土柱凍脹量隨凍結時間不斷增大，最終凍脹量約為5 mm，數值模擬結果與實測結果趨勢變化一致，驗證了該模型的有效性。