多年凍土區(qū)路基填筑短期內(nèi)凍土溫度場數(shù)值分析

王蘊嘉，郭惠芹，葉陽升，王仲錦，王鵬程，張千里

（1.中國鐵道科學研究院集團有限公司鐵道建筑研究所，北京 100081；2.中國鐵道科學研究院研究生部，北京 100081；3.中國鐵道科學研究院集團有限公司，北京 100081）

我國青藏高原的平均海拔超過4 500 m，其中近70%的面積被多年凍土覆蓋［1-2］。多年凍土是一種對溫度變化非常敏感的特殊土體，其物理力學性質(zhì)與土體溫度密切相關(guān)。在多年凍土上進行路基修建等活動，勢必會打破凍土與環(huán)境間原有的能量動態(tài)平衡，改變凍土場的熱狀態(tài)，進而影響路基等構(gòu)筑物的穩(wěn)定性［3-6］。青藏鐵路穿越多年凍土區(qū)路段長達632 km，在建的川藏鐵路同樣面臨多年凍土問題?；诖耍_展多年凍土區(qū)路基填筑對凍土溫度場影響的相關(guān)研究對保證兩路的順利修建和安全運營具有重要意義。

現(xiàn)場監(jiān)測和試驗是研究路基填筑對凍土溫度場影響的重要方法。王小軍等［7］分析了11 個青藏鐵路路基斷面7 年監(jiān)測數(shù)據(jù)，發(fā)現(xiàn)路基中心人為上限均隨時間發(fā)展有所上升，有利于路基穩(wěn)定；牛富俊等［8］調(diào)研了青藏鐵路沿線55 個斷面的監(jiān)測結(jié)果，發(fā)現(xiàn)部分高溫凍土區(qū)普通路基存在凍土上限下降、路基熱穩(wěn)定性差的情況，且工程施工對地表的熱擾動、水熱侵蝕易引發(fā)熱融性災害，產(chǎn)生路基沉陷等病害；穆彥虎等［9］分析了2005 年—2010 年間青藏鐵路普通路基斷面的地溫監(jiān)測數(shù)據(jù)，發(fā)現(xiàn)多年凍土區(qū)路基下凍土熱穩(wěn)定性與多年凍土平均地溫密切相關(guān)，高溫凍土區(qū)路基下部凍土熱穩(wěn)定性較差；張明禮等［10］研究發(fā)現(xiàn)，路基運營初期陽坡凍土上限下降，水的滲流不僅加速了路基的升溫和變形過程，還加劇了路基溫度場和水分場的非均勻性。此外，趙相卿等學者［11-12］也進行了相關(guān)研究。然而，現(xiàn)有研究大多關(guān)注的是路基填筑對多年凍土長期熱穩(wěn)定性的影響，而對路基填筑后短期內(nèi)的溫度場變化規(guī)律研究較少。且對短期內(nèi)的凍結(jié)鋒面下凹扭曲的監(jiān)測存在一定難度。實際工程中，路基填筑熱效應在短期內(nèi)使得路基覆蓋范圍凍結(jié)鋒面下凹扭曲，且在水分匯入充足的情況下還將引起持續(xù)性的下凹扭曲，產(chǎn)生路基熱融沉降、路基開裂等病害，有必要對其進行深入研究。

近年來，大量學者采用數(shù)值模擬研究凍土工程相關(guān)問題，合理有效的凍土水熱耦合模型則是數(shù)值研究的關(guān)鍵。Harlan［13］基于水分遷移理論提出了水熱耦合的概念，首次建立了滲透系數(shù)與含冰量相關(guān)的水熱耦合模型；徐學祖等［14］基于大量試驗研究，給出了未凍水含量的經(jīng)驗公式；Liu 和Tian 等學者［15-16］以土體的狀態(tài)焓為變量模擬分析相變潛熱的影響；宋二祥等［17］基于Philip 等提出的非等溫水熱氣耦合運移模型，提出綜合考慮水力梯度和溫度梯度影響的水氣運移控制方程。此外，白清波等［18-19］也對凍土水熱耦合模型進行修正和改進，并將其應用于數(shù)值模擬研究中。

基于凍土水熱耦合計算控制方程，采用有限元軟件COMSOL 提供的系數(shù)型偏微分方程模塊（PDE）建立凍土水熱耦合分析模型，進行路基填筑短期內(nèi)凍土溫度場的變化規(guī)律研究。

1 凍土水熱耦合計算控制方程

1）水分場控制方程

凍結(jié)和融化土體中均含有一定量的未凍水。根據(jù)Philip 等［20］提出的非等溫水熱氣耦合運移模型，水分遷移主要由水力梯度和溫度梯度導致。參考文獻［17］提出的水氣運移控制方程，不考慮氣態(tài)水的遷移，非飽和土中液態(tài)水遷移的控制方程可表示為

式中：θu和θi分別為未凍水和冰的體積分數(shù)；t 為時間，s；ρi為冰的密度，取900 kg·m-3；ρw為水的密度，取1 000 kg·m-3；Kh為等溫條件下水力梯度導致的水力傳導系數(shù)，m·s-1；KT為非等溫條件下溫度梯度導致的水力傳遞系數(shù)，m2·K-1·s-1；h 為基質(zhì)吸力對應水頭，m；hg為重力作用引起的水頭，即位置水頭，m；T 為絕對溫度，K；?為微分算子。

其中，

m=1-1/n

式中：α，n 和m 均為模型參數(shù)，需通過試驗確定，α 單位為m-1；θs和θr分別為飽和體積含水率和殘余體積含水率。

考慮冰水混合區(qū)冰對液態(tài)水遷移的阻礙作用，Kh可表示為

式中：ks為土的滲透系數(shù)。

非等溫條件下溫度梯度導致的水力傳遞系數(shù)KT的計算公式［21］為

2）溫度場控制方程

忽略土體凍結(jié)和融化過程中水分蒸發(fā)耗熱等的影響，將冰水相變潛熱視為熱源，根據(jù)傅里葉定律，非飽和土瞬態(tài)溫度場的控制方程可表示為

式中：C 為土的體積熱容，J·（m3·K）-1；λ 為土的導熱系數(shù)，W·m-1·K-1；L為冰水相變潛熱，取334.5 kJ·kg-1。

凍土的體積熱容和導熱系數(shù)均由土骨架、未凍水、冰和孔隙氣體共同決定。綜合考慮計算效率和收斂性，本文僅針對土的凍融狀態(tài)給定不同的熱力學參數(shù)，即

式中：Cf和Cuf分別為凍土和融土的體積熱容，J·m-3·K-1；λf和λuf分別為凍土和融土的導熱系數(shù)，W·m-1·K-1；Tf為土的凍結(jié)溫度，K。

3）未凍水含量

水分場控制方程和溫度場控制方程中共包含3 個未知量，分別為未凍水體積分數(shù)、冰體積分數(shù)和溫度。若想求解水熱耦合問題，仍需引入1個關(guān)于3個未知量的聯(lián)系方程。

參考文獻［19］未凍水含量可由下式計算。

式中：B為經(jīng)驗參數(shù)。

2 凍土水熱耦合計算控制方程驗證

本文采用有限元軟件COMSOL，模擬文獻［14］中封閉條件下的土柱融化試驗，進行驗證凍土水熱耦合計算控制方程的有效性。

文獻［14］的融化試驗采用蘭州粉土制成的土柱直徑10 cm、高度10 cm，初始體積含水率ω 為22.3%，初始溫度為-1 ℃，試驗過程始終保持土柱底部溫度為-1 ℃，頂部溫度為1 ℃。

采用有限元軟件COMSOL 建立高度為10 cm的一維土柱模型，土體采用二次拉格朗日單元模擬，水分場、溫度場的初始值和邊界條件與文獻［14］的試驗相同。參考文獻［17-18］，模擬中粉土參數(shù)取值見表1。文獻［14］的試驗和數(shù)值模擬得到的試樣含水率分布、試樣溫度分布如圖1和圖2所示。遷移，所以在凍結(jié)溫度附近水分大量聚集，形成含水率峰值；土柱溫度基本成線性分布，符合實際。

表1 融化試驗模型參數(shù)取值

圖1 融化試驗中的體積含水率分布

圖2 融化試驗中試樣的溫度分布

可見，采用本文凍土水熱耦合計算控制方程可較準確地模擬凍土中水分遷移和溫度場變化，可采用該方法對多年凍土區(qū)路基填筑問題進行模擬研究。

3 模型建立

參考青藏鐵路北麓河試驗段路基斷面［19-20］，采用有限元軟件COMSOL 建立路基模型。地基寬50 m，深20 m；路基寬8 m，高4 m，坡比固定為1∶1.5。為提高計算效率，取對稱模型的一半進行計算，如圖3 所示。路基和地基均采用三角形單元模擬，主要計算參數(shù)見表2。

表2 路基模型主要參數(shù)

圖3 模型示意圖

由圖1 和圖2 可知：模擬結(jié)果和試驗實測值基本吻合；含水率在土柱上半部分基本恒定且小于初始值，這是上部凍土融化后液態(tài)水在水力梯度和溫度梯度共同作用下向下遷移導致；含水率在土柱下半部分隨溫度降低先增大后減小，出現(xiàn)含水率峰值，這是因為凍結(jié)溫度附近水力梯度大、水分遷移速度快，而冰的阻滯作用使得液態(tài)水難以繼續(xù)向下

參考文獻［8，22-23］中青藏鐵路北麓河試驗段相關(guān)監(jiān)測數(shù)據(jù)確定模型的溫度邊界。模型溫度場左、右邊界為絕熱邊界，底部采用狄利克雷邊界，取-0.8 ℃，頂部溫度Ttop隨時間周期變化，即

由于高原高寒地區(qū)一般在暖季施工，即每年的4 月至10 月，因此式（9）中的t=0 s 時對應時間為每年的4月15日。模型水分場左、右邊界為零通量邊界，底部、頂部均為狄利克雷邊界，底部未凍水體積分數(shù)為0.12，頂部考慮蒸發(fā)作用較為干燥，未凍水體積分數(shù)為0.05。

地基土的水分場和溫度場初始值通過計算求得。在沒有路基條件下，給定地基土初始溫度和未凍水體積分數(shù)分別為-0.8 ℃和0.12，然后采用前文所述邊界條件進行計算，求解得到50 a后地基土的溫度場和水分場初始值。

經(jīng)過運輸、壓實等過程的填料溫度往往高于地表平均溫度，偏安全考慮，路基填筑完填料溫度比地表溫度高10 ℃。填料未凍水體積分數(shù)為0.12。本文重點關(guān)注短期內(nèi)路基填筑熱效應，計算時長取為90 d。

4 路基填筑短期內(nèi)凍土溫度場

路基的填筑會在短期內(nèi)對凍土溫度場產(chǎn)生較大影響，影響因素主要來源于初始條件、路基尺度和填料熱力學性質(zhì)3方面。初始條件包括路基填筑時間和填料與地表溫差。對于地基而言，4 月至7 月地表溫度持續(xù)升高，地基土從地表吸入熱量；7 月至10 月地表溫度逐漸減低，地基土向地表輸出熱量。因此，分別選取4 月、7 月和10 月3 個關(guān)鍵時間節(jié)點進行計算，研究填筑時間的影響。路基尺度包括路基高度和寬度。填料熱力學性質(zhì)包括體積熱容和導熱系數(shù)。

4.1 填筑時間的影響

圖4 和圖5 分別為4 月、7 月和10 月進行路基填筑后30 和60 d 的溫度場，圖中紅色粗實線代表凍結(jié)鋒面，即0 ℃刻度線位置。

由圖4 和圖5 可知：填筑后30 d 時，最高溫度出現(xiàn)在路基中心位置，說明此時路基填料的熱量仍在持續(xù)向大氣和地基傳遞；填筑后60 d 時，4 月填筑路基溫度場分布從地表到地基逐漸降低，與天然地基的溫度分布規(guī)律一致，說明路基填料的熱量已經(jīng)基本耗散完畢，而7月和10月填筑路基高溫區(qū)域均出現(xiàn)在路基中心位置且顯著高于天然地基，說明填料的熱量仍在持續(xù)向大氣和地基傳遞。

圖4 填筑后30 d的溫度場

圖5 填筑后60 d的溫度場

冰的阻滯作用使得凍結(jié)土體滲透系數(shù)顯著低于融化土體，凍土融化產(chǎn)生的水分以及從地表滲入地基的水分往往聚集在凍結(jié)鋒面附近難以下滲，并在重力的作用下沿凍結(jié)鋒面向下凹處流動、匯集。從地表滲入的水分溫度較高，帶有一定熱量，在下凹處匯集會使得該處凍土受到熱侵蝕繼續(xù)融化、凍結(jié)鋒面下降、下凹程度增大，而這促使周圍水分繼續(xù)向凍土鋒面下凹處匯集，使得該處凍結(jié)鋒面下凹扭曲程度持續(xù)增大。

實際工程中，路基填筑的熱效應使得路基覆蓋范圍凍結(jié)鋒面下凹扭曲，即便扭曲程度較小、與天然凍結(jié)鋒面差值僅幾個厘米，在水分匯入充足的情況下也可能引起持續(xù)性的下凹扭曲，產(chǎn)生路基熱融沉降、路基開裂等病害。因此，相較于凍結(jié)鋒面的絕對高程，更關(guān)心路基填筑后的凍結(jié)鋒面下凹扭曲程度，即路基覆蓋范圍凍結(jié)鋒面與天然地基凍結(jié)鋒面的高程差值。

圖6 為路基覆蓋范圍凍結(jié)鋒面高程差隨時間發(fā)展曲線，正值表示該處凍結(jié)鋒面位置高于天然場地，負值表示低于天然場地，即該處凍結(jié)鋒面發(fā)生下凹扭曲。

圖6 路基覆蓋范圍凍結(jié)鋒面高程差隨時間發(fā)展曲線

由圖6（a）可知：4 月填筑時，路基中心凍結(jié)鋒面先下凹扭曲，最大下凹深度20 cm，然后隨時間增大下凹程度逐漸減小，最后變成上凸扭曲；坡腳位置填筑后持續(xù)上凸。從工程角度出發(fā)，坡腳凍結(jié)鋒面上凸扭曲可以阻礙路基外側(cè)水的匯入，對工程而言是有利的。但同時坡腳上凸可能會使得路基外側(cè)水分以及沿坡面流下的水分在坡腳外側(cè)聚集，引發(fā)坡腳處的水熱侵蝕。

由圖6（b）可知：7 月填筑時，路基覆蓋范圍凍土場下凹扭曲程度持續(xù)增大，路基中心凍結(jié)鋒面高程差從0 增至25 cm；路基坡腳僅在填筑后90 d才觀察到輕微上凸扭曲。夏季降雨較多，地表水豐富且溫度高、熱量大，此時凍土場發(fā)生下凹扭曲最容易引發(fā)水熱侵蝕，于工程而言是最危險的。

由圖6（c）可知：10 月填筑時，路基覆蓋范圍凍土場下凹扭曲程度持續(xù)增大，路基中心凍結(jié)鋒面高程差從0 增至38 cm；坡腳并未觀察到上凸扭曲。雖然10 月填筑后相同時間內(nèi)的凍土場下凹扭曲程度最大，但考慮到填筑后地表溫度很快降到0 ℃以下，地表水較少且溫度低，不易引起水熱侵蝕。

圖7 為4 月、7 月和10 月填筑后路基中心與天然場地凍結(jié)鋒面高程差隨時間發(fā)展曲線。由圖7可知：4 月填筑時路基中心凍結(jié)鋒面高程差發(fā)展規(guī)律與7 月、10 月填筑不同，4 月填筑后路基中心凍結(jié)鋒面先下降后上升，而7月、10月填筑后路基中心凍結(jié)鋒面持續(xù)下降，即路基覆蓋范圍下凹扭曲程度持續(xù)增大；7 月填筑后凍結(jié)鋒面高程差絕對值與時間近似成對數(shù)關(guān)系，而10 月填筑后凍結(jié)鋒面高程差隨時間增加而急劇增大，直至路基覆蓋范圍形成封閉的融化核。后續(xù)研究將以7 月填筑為基礎展開。

圖7 路基中心凍結(jié)鋒面高程差隨時間發(fā)展曲線

4.2 路基填料與地表溫差的影響

圖8 和圖9 分別為填料溫度與地表溫度差不同時，填筑后90 d路基覆蓋范圍凍結(jié)鋒面下凹扭曲曲線和路基中心與天然場地凍結(jié)鋒面高程差。

圖8 路基覆蓋范圍凍結(jié)鋒面高程差曲線

圖9 路基中心凍結(jié)鋒面高程差隨路基填料與地表溫度變化曲線

由圖8 可知：路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨填料與地表溫差增大而增大；當溫差為4 ℃時，在坡腳位置還可觀察到明顯的凍結(jié)鋒面上凸扭曲，可在一定程度上阻礙周圍水分向路基中心底部匯聚。

任一超過材料或結(jié)構(gòu)疲勞極限的應力的循環(huán)作用均會對材料或結(jié)構(gòu)造成不同程度的損傷，而造成的損傷是累積的，當損傷累積達到材料或結(jié)構(gòu)疲勞失效的臨界值時材料或結(jié)構(gòu)就會出現(xiàn)疲勞損壞。

由圖9 可見：路基中心與天然場地凍結(jié)鋒面高程差隨溫差增大而增大。

4.3 路基尺度的影響

圖10 和圖11 分別為路基高度不同時，填筑后90 d路基覆蓋范圍凍結(jié)鋒面下凹扭曲曲線和路基中心與天然場地凍結(jié)鋒面高程差。

圖10 路基覆蓋范圍凍結(jié)鋒面高程差曲線

圖11 路基中心凍結(jié)鋒面高程差隨路基高度變化曲線

由圖10 可知：路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨路基高度增大而增大。

由圖11 可見：路基中心與天然場地凍結(jié)鋒面高程差隨路基高度增大而增大。

圖12 和圖13 為路基寬度不同時，填筑后90 d路基覆蓋范圍凍結(jié)鋒面下凹扭曲曲線和路基中心與天然場地凍結(jié)鋒面高程差。

由圖12 可知：路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨路基寬度增大而增大。

圖12 路基覆蓋范圍凍結(jié)鋒面高程差曲線

由圖13 可見：路基中心與天然場地凍結(jié)鋒面高程差隨路基寬度增大而增大。

圖13 路基中心凍結(jié)鋒面高程差隨路基寬度變化曲線

4.4 填料熱力學性質(zhì)的影響

圖14 和圖15 分別為填料導熱系數(shù)（凍結(jié)狀態(tài)）不同時，填筑后90 d路基覆蓋范圍凍結(jié)鋒面下凹扭曲曲線和路基中心與天然場地凍結(jié)鋒面高程差。

圖14 路基覆蓋范圍凍結(jié)鋒面高程差曲線

圖15 路基中心凍結(jié)鋒面高程差隨填料導熱系數(shù)變化曲線

由圖14 可知：路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨填料導熱系數(shù)增大而減小；當填料導熱系數(shù)為2.0 W·m-1·℃-1時，在坡腳位置還可觀察到明顯的凍結(jié)鋒面上凸扭曲，可在一定程度上阻礙周圍水分向路基中心底部匯聚。

由圖15 可見：路基中心與天然場地凍結(jié)鋒面高程差隨填料導熱系數(shù)增大而減小。

圖16和圖17分別為填料體積熱容（凍結(jié)狀態(tài)）不同時，填筑后90 d路基覆蓋范圍凍結(jié)鋒面下凹扭曲曲線和路基中心與天然場地凍結(jié)鋒面高程差。

圖16 路基覆蓋范圍凍結(jié)鋒面高程差曲線

由圖16 可知：路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨填料體積熱容增大而增大。

由圖17 可見：路基中心與天然場地凍結(jié)鋒面高程差隨填料體積熱容增大而增大。

5 路基填筑短期內(nèi)凍結(jié)鋒面計算方法

凍結(jié)鋒面下凹扭曲曲線可用倒梯形近似描述，如圖18所示。

圖18 凍結(jié)鋒面下凹扭曲示意圖

該曲線可分為AB，BC和CD 3段：AB段為路基中心部分，該區(qū)域與天然場地凍結(jié)鋒面高程差最大，且差值基本恒定，不隨水平位置的變化而改變；BC段為路基肩部到路基坡腳的部分，該區(qū)域與天然場地凍結(jié)鋒面高程差迅速降低，近似線性變化；CD段對應路基坡腳及相鄰區(qū)域，該區(qū)域與天然場地凍結(jié)鋒面高程差基本為0。因此凍結(jié)鋒面高程差Δh 可表示為

式中：Δhmax為路基中心的凍結(jié)鋒面高程差，即圖18 中的OA；W1為路基中心的凍結(jié)鋒面高程差穩(wěn)定段寬度，即圖18 中的AB；W2為凍結(jié)鋒面變化寬度，即圖18 中的OC；W 為路基中心與坡腳距離，即圖18中的OD。

以上參數(shù)的具體計算如下。

（1）Δhmax與時間近似成對數(shù)關(guān)系，隨填料與地表溫差、路基高度、路基寬度和填料體積熱容的增大而增大，隨填料導熱系數(shù)增大而減小。根據(jù)圖7、圖8、圖9、圖10、圖11、圖13 和圖15 中的計算結(jié)果，擬合得到Δhmax為

式中：tf為填筑后時長，d；ΔT 為填筑時填料與地表溫差，℃；h 為路基高度，m；W路基為路基寬度，m；Cf為填料凍結(jié)時的體積熱容，J·m-3·K-1；λf為填料凍結(jié)時的導熱系數(shù)，W·m-1·K-1。

（2）W 由路基高度和寬度2個參數(shù)確定，計算公式如下。

（3）W1由路基高度和寬度2 個參數(shù)確定，計算公式如下。

（4）W2隨路基高度和寬度的增大而增大，隨填料導熱系數(shù)的增大而減小，根據(jù)圖14 中的計算結(jié)果，擬合得到計算式如下，即

6 結(jié) 論

（1）路基填筑后短期內(nèi)路基中心下方凍結(jié)鋒面發(fā)生下凹扭曲，在水分匯入充足情況下將引發(fā)水熱侵蝕，導致凍結(jié)鋒面持續(xù)性的下凹扭曲，進而可能產(chǎn)生路基熱融沉降、路基開裂等病害，工程中需做好防水工作。

（2）與4 月填筑和10 月填筑相比，7 月填筑時最大程度阻礙了土體熱量向大氣中輸出，路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度最大、凍土溫度場恢復最慢，屬于最不利工況。

（3）路基覆蓋范圍凍結(jié)鋒面下凹扭曲程度隨填料與地表溫差、路基高度、路基寬度、填料體積熱容的增大而增大，隨填料導熱系數(shù)增大而減小，因此實際工程中應盡量選用體積熱容小、初始溫度低，且導熱系數(shù)大的填料。

（4）提出路基覆蓋范圍凍結(jié)鋒面下凹扭曲計算方法，并給出關(guān)鍵參數(shù)確定方法。采用該方法可得到填筑完成后不同時刻、不同工況下的凍結(jié)鋒面下凹扭曲曲線。