寒區(qū)土體一維水-熱-力耦合模型與數(shù)值分析

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

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

引言

凍土在我國國土面積中占據(jù)著極其廣闊的一部分，其中多年凍土約占1/4。寒區(qū)凍土的凍脹融沉特性使結(jié)構(gòu)物產(chǎn)生不均勻沉降、結(jié)構(gòu)斷裂等現(xiàn)象，工程災(zāi)害頻頻發(fā)生，對上覆結(jié)構(gòu)物的穩(wěn)定帶來了極大的挑戰(zhàn)。凍土凍融作用的機(jī)理是當(dāng)前凍土研究中最受關(guān)注的問題之一［1-5］。

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

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

1 水-熱-力三場耦合數(shù)學(xué)模型

1.1 基本假設(shè)

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

（1）在土體中，水分遷移主要以液態(tài)形式進(jìn)行遷移，在模型中忽略氣態(tài)水的遷移；

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

（3）土體各向同性，均勻連續(xù)；

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

1.2 水分場控制方程

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

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

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

則單位時間內(nèi)土體微單元的水分增量為：

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

聯(lián)立方程可得：

式中ρi為冰的密度。

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

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

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

1.3 溫度場控制方程

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

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

熱通量q又為導(dǎo)熱系數(shù)λ與溫度T的函數(shù)：

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

聯(lián)立方程建立溫度場控制方程：

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

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

文中水分場與溫度場之間屬強(qiáng)耦合，彼此相互影響。方程（15）以土體凍結(jié)溫度Tf為分界點可分為2個階段：當(dāng)溫度超過土體的凍結(jié)溫度時，土體內(nèi)不會發(fā)生冰水相變過程，此時溫度變化與土體的熱容以及導(dǎo)熱系數(shù)有關(guān)；當(dāng)溫度低于土體的凍結(jié)溫度時，土體內(nèi)水開始相變成冰。一方面，相應(yīng)位置含水量減少，基質(zhì)勢發(fā)揮作用，土體內(nèi)水分開始遷移，由液態(tài)水含量高的位置遷移至液態(tài)水含量低的位置。另一方面，冰水相變所產(chǎn)生的潛熱也開始影響溫度的變化速率，水分場與溫度場之間的耦合開始發(fā)揮作用。

1.4 應(yīng)力場控制方程

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

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

靜力平衡方程：

幾何方程：

物理方程：

在凍結(jié)過程中，凍結(jié)區(qū)內(nèi)的水分產(chǎn)生了相變，水相變成冰體積相應(yīng)增大9%。與此同時，在基質(zhì)勢的作用下，未凍結(jié)區(qū)的液態(tài)水向凍結(jié)區(qū)移動，則土體的應(yīng)變可以用遷移的液態(tài)水含量來表示，即溫度作用下由于水分遷移產(chǎn)生的應(yīng)變εvf為：

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

需要說明的是，土體凍融是一個非常復(fù)雜的物理過程，溫度場以及水分場變化影響應(yīng)力的重分布。反之，應(yīng)力重分布又影響著土體溫度以及體積含水量的變化［16］。在溫度場作用下，凍結(jié)區(qū)的水分發(fā)生相變導(dǎo)致凍結(jié)區(qū)液態(tài)水含量降低，進(jìn)而在基質(zhì)勢作用下未凍結(jié)區(qū)液態(tài)水遷移至凍結(jié)區(qū)，其填充土體內(nèi)孔隙并產(chǎn)生進(jìn)一步相變，水相變成冰體積增大，土體內(nèi)部結(jié)構(gòu)產(chǎn)生應(yīng)力。在應(yīng)力作用下土體發(fā)生應(yīng)變，進(jìn)而造成土體的凍脹。

2 基于Comsol數(shù)值模型構(gòu)建

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

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

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

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

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

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

溫度場控制方程為：

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

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

圖4 Comsol內(nèi)置假設(shè)方程形式Fig.4 Equation form in Comsol

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

表2 數(shù)值模型參數(shù)表Table 2 Table of numerical model parameters

3 結(jié)果與討論

3.1 溫度場

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

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

圖6 距離土柱頂部5 cm處溫度隨凍結(jié)時間變化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

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

3.2 水分場

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

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

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

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

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

3.3 凍脹位移

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

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

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

4 結(jié)論

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

文中基于質(zhì)量守恒定律以及能量守恒定律對凍土水分場以及溫度場控制方程進(jìn)行系統(tǒng)的推導(dǎo)，并結(jié)合水熱條件下凍土內(nèi)部應(yīng)力變化建立水熱力三場耦合數(shù)學(xué)模型，通過建立一維非飽和土柱數(shù)值模型可得結(jié)論：

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

（2）土柱內(nèi)水分在凍結(jié)過程中發(fā)生遷移，在凍結(jié)面處形成明顯的S型曲線，未凍結(jié)區(qū)水分不斷向凍結(jié)區(qū)遷移；

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