水-熱耦合作用下嚴寒地區高速鐵路隧道溫度場分布規律

趙希望, 張 恒

(1.廣東省隧道結構智能監控與維護企業重點實驗室， 廣州 511440； 2.中鐵隧道局集團有限公司勘察設計研究院， 廣州 511440； 3.華南理工大學土木與交通學院， 廣州 510641)

寒冷地區隧道的修筑打破了圍巖的原始熱平衡，圍巖中存在的水分會在寒季凍結成冰[1]。由于液相轉為固相體積增大，進而使隧道襯砌背后圍巖產生凍脹力；在暖季，凍結的冰融化成水，在凍結時產生的位移不能完全恢復，造成的空隙又會增加滲漏，嚴重時可能發生融陷，這種襯砌背后圍巖的反復凍融作用使襯砌防水層發生破損，進而產生更大的凍害，造成襯砌結構變形，使隧道運營安全受到威脅[2-4]。

眾多專家和學者對新疆、青海等嚴寒地區的隧道進行了研究。賴遠明等[5]運用伽遼金法對青海大坂山隧道溫度場和滲流場耦合非線性問題進行了分析和研究，通過無量綱分析法和攝動技術求出寒區圓形隧道溫度場的解析解。高焱等[6]對寒區隧道洞內空氣溫度及圍巖溫度場進行了分析，給出了保溫層的鋪設范圍。陳建勛等[7]首先結合寒區公路隧道防凍保溫層的設計，對保溫材料現場的保溫性能進行了測試，確定了最佳的防凍保溫材料組合和施工，然后運用正弦函數回歸法對隧道拱頂、拱腰、邊墻和路面四個部位實測的溫度數據進行了分析。張學富等[8]根據青藏鐵路昆侖山隧道內實測大氣溫度和圍巖溫度，考慮冰-水相變影響，應用伽遼金法推導了有限元公式，對昆侖山隧道施工期間融化的圍巖進行了回凍預測研究，分析了隧道鋪設保溫材料和無保溫材料的溫度場，得出保溫材料對多年凍土區隧道融化回凍過程所起的作用。郭春香等[9]以昆侖山隧道為例，利用有限元法分別計算了考慮混凝土水化熱和不考慮混凝土水化熱的隧道圍巖融化回凍過程，分析了混凝土水化熱對圍巖融化、回凍過程的影響。楊旭等[10]基于圍巖溫度場現場測試的結果，同時考慮水文地質條件、混凝土水化熱、大氣溫度和地溫等因素，比較了季節凍土區隧道施加保溫層和未施加保溫層的凍融循環圈的差異。周元輔等[11-12]推導了多年凍土隧道隔熱層的數學優化模型，分析了隔熱層用于多年凍土隧道的合理性和經濟性。從以上對寒區隧道的研究中可以看出，防凍保溫層是寒區隧道防治凍害的有效措施，確定保溫層鋪設的合理性和經濟性對寒區隧道工程有著重要的意義。但是，以上分析中或是以多年凍土區隧道施工期間造成的圍巖融化為對象，研究襯砌背后圍巖的回凍過程，并不適用于初始融土地區隧道修筑后，由于圍巖熱平衡破壞造成的襯砌背后圍巖凍結的問題，如文獻[8-9]；或是在隧道溫度場的分析中只考慮了襯砌、圍巖熱傳導，沒有考慮水分遷移或冰水相變，如文獻[10-12]。

針對融土地區隧道修筑后，襯砌背后圍巖凍結的問題，以吉圖琿客運專線榆樹川隧道為依托，根據瞬態溫度場問題的能量平衡方程和水分遷移方程，考慮水分遷移、冰-水相變等因素，建立了寒區隧道圍巖水-熱耦合模型，并結合現場實測大氣溫度，利用有限元法對嚴寒地區隧道溫度場進行計算，比較有、無保溫層對隧道溫度的影響，確定隧道凍害發生的最不利位置和時間，分析保溫材料對隧道凍融圈的影響，給出保溫層厚度和導熱系數的合理值。

1 水-熱耦合作用下傳熱數學模型

1.1 溫度場方程

土體在凍結時水蒸氣蒸發耗熱很少，考慮水分遷移和冰水相變，土體內發生的熱傳導[13]可以表示如下。

融化區:

(1)

凍結區:

(2)

式中：f、u為圍巖凍結和融化狀態；λ為熱傳導系數，W/(m·℃)；T為溫度，℃；t為時間；c為圍巖的比熱容，J/(kg·℃)；ρ和ρi為分別為圍巖的密度和冰的密度，kg/m3；L為單位質量冰水相變潛熱，kJ/kg；θi為體積含冰量。

在凍結鋒面位置處s(t)，應滿足連續性條件

Tf[s(t),t]=Tu[s(t),t]=Tm

(3)

式(3)中：Tm為土體的凍結溫度，℃。

1.2 水分場方程

根據文獻[14-15]，土體內的空氣和水蒸氣在孔隙內的流動帶動的水分遷移量相當微弱，因此，忽略空氣和水蒸氣遷移對水分的遷移的影響，凍結和融化過程中圍巖內水分遷移滿足Darcy定律，可表達為

(4)

式(4)中：K為導水系數，m/s；φ為基質勢，Pa；θu為未凍水體積含量；ρw為水的密度，kg/m3。

根據徐斅祖等[16]的實驗結果，未凍水含量θu與負溫始終保持動態平衡的關系，即

θu=a1|T|-b1

(5)

式(5)中：a1和b1為與土質因素有關的經驗常數。

1.3 溫度場方程和水分場方程的耦合

土體中的微分水容量C(m-1)是等溫條件下基質勢與體積含水量θu之間的定量關系[16]，表達式為

(6)

土壤的水分擴散系數D可以表示為土體微分水容量C與導水系數K的關系[16]，即

(7)

Taylor等[17]認為未凍結土的土-水特征曲線也適用于凍結土的情況。因此，根據凍土中未凍水含量與溫度的關系以及土-水特征曲線可確定不同溫度的未凍水基質勢，表達式為

(8)

綜合式(2)、式(4)、式(6)～式(8)，土體凍結區的熱傳導方程可以表示為

(9)

用等效熱傳導系數λe和等效比熱容Ce代替圍巖融化和凍結狀態的熱傳導系數λu、λf和比熱容cu、cf，則考慮冰-水相變和水分遷移的圍巖熱傳導方程可以表示為

(10)

式(1)～式(10)組成了隧道圍巖凍結過程中的水-熱耦合模型，采用COMSOL Multiphysics 軟件的PDE (partial differential equations)模塊進行求解。以T和θ為因變量，得到兩個廣義型PDE標準形式為

(11)

式(11)中：da為質量系數；U為微分方程因變量，U={T,θ}；Γ為通量向量；f為源項。

將節點溫度T和含水量θ作為基本未知量，運用隱式歐拉向后差分(BDF)法對控制方程中的時間項進行離散，采用非線性迭代修正阻尼牛頓法求解方程。利用有限差分時間步長的方法，對該微分方程進行時間差分，并通過連續線性循環解法求解。

1.4 模型驗證

為驗證提出的水-熱耦合模型的正確性，采用胡和平等[18]對于張掖壤土的凍結試驗和結果進行數值模擬。胡和平等采用張掖壤土，在試驗室內進行了一組單向凍結試驗。試樣的高度為13.68 cm，直徑為11.36 cm。試樣的初始溫度如表1所示，試樣的初始體積含水量為22.08%。凍結實驗時，試樣的側面為熱絕緣條件，頂部與底部溫度按照如圖1所示變化，試驗在凍結2 830 min后結束，并將試樣取出，然后將試樣切片，在烤箱中烘干得到總含水量。

圖1 試樣頂部和底部溫度變化Fig.1 Temperature variations at the top and bottom of the sample

圖2 計算模型Fig.2 Calculation model

表1 試樣初始溫度分布

計算時，取軸對稱計算模型(圖2)，邊界條件和初始條件如下：

溫度邊界條件：AC為對稱邊界；BD為絕緣邊界；初始溫度場沿深度方向按照表1分布。AB和CD邊界溫度隨時間變化情況如圖1所示。

滲流邊界條件：AC為對稱邊界；AB、BD和CD都是不透水邊界，即試樣處于封閉環境中。

試驗所用土為張掖壤土，熱物理參數如表2所示。

表2 張掖壤土熱物理參數

在土體凍結過程中，土體孔隙中生成的冰晶體會阻礙孔隙內水分的遷移孔道，從而導致水分擴散系數D和滲透系數K減小。受土體孔隙內冰晶體的影響，土體在凍結過程中的分水擴散系數D和滲透系數K[19]可以表示為

(12)

(13)

式中：I為阻抗系數；a2、b2、a3和b3分別為物理試驗常數，其值如表3所示。

表3 張掖壤土物理參數

圖3所示為試樣在不同凍結時刻(34、528、2 830 min)，不同深度處溫度實驗值與計算值的比較。可以看出，各凍結時刻試樣內溫度的計算值和實驗值能夠很好地相吻合，這說明該水-熱模型能夠很好地模擬土體凍結過程中的溫度分布規律。同時，可以看出：凍結初始階段，試樣處于正溫狀態，但隨著凍結的進行，試樣發生自上而下的凍結，凍結鋒面不斷下移。凍結初期溫度傳遞速率較快，隨凍結時間的增加溫度傳遞速率逐漸減小，最終溫度趨于穩定。

圖3 不同時刻土樣溫度計算值與試驗值的比較Fig.3 Comparison between the experimental and simulated temperatures of soil sample at different time

圖4為實驗結束時試樣內含水量的計算值和實驗值的比較，含水量的實驗值和計算值在形狀和分布特征上都是一致的。試驗開始時，試樣的初始含水量為22.08%，但試驗結束時，試樣內的最大含水量的實驗值為29.2%，計算值為32.2%；最小值含水量的實驗值為18.2%，計算值為18.4%，計算值與實驗值非常接近，并且最大含水量位置的計算值和實驗值也相當吻合。這說明水-熱模型可以很好地模擬土體凍結過程中的水分遷移現象。

圖5為在不同凍結時刻(100、528、2 830 min)，試樣不同深度處的含冰量分布。可以看出，在凍結過程中試樣發生自上而下的凍結，試樣含冰量隨時間的變化也表現出凍結初期主要分布在土體淺層，隨著凍結的進行，冰向土體內部不斷增長，在凍結鋒面位置處含冰量達到最大。由于在實驗中土體含冰量的變化無法進行測量，在實驗結束后測得的都是液態水的含量，所以對土柱在凍結過程中冰含量的分布和變化的分析是一個難題。水-熱耦合模型可以實現對土體在凍結過程中冰含量變化的監測，為分析土體內冰的分布和含量提供了方便。

圖4 實驗結束時土樣含水量的計算值和實驗值的比較Fig.4 Comparison between the experimental and simulated water content of soil sample at the end of experiment

圖5 不同時刻土體含冰量的計算值Fig.5 The calculated ice content of soil sample at different time

2 數值模擬

2.1 工程概況

新建吉圖琿客運專線位于吉林省東部吉林市和延邊朝鮮族自治州境內，線路全長360 km，設計時速250 km/h，其中，榆樹川隧道進口位于延吉市安圖縣境內附近，隧道全長2 210 m，最大埋深約160 m，如圖6所示。根據延吉市地區氣象資料和水文地質條件，隧址年平均氣溫5.32 ℃，極端最高氣溫37.7 ℃，極端最低氣溫-29.2 ℃，最冷月平均氣溫-16.5 ℃，最大凍結深度為1.68 m。隧道進口處大氣月平均溫度如圖7所示。

圖6 榆樹川隧道照片Fig.6 Photo of the Yushuchuan tunnel

圖7 隧道外大氣溫度Fig.7 Air temperature outside the tunnel

根據水文地質條件，榆樹川隧道通過地層為強風化、全風化砂巖。隧道開挖后，地層之間的熱平衡被打破，隧道的進口和出口處因受到大氣溫度的作用，在冷季圍巖中的水分發生凍結，從而使隧道襯砌發生凍脹破壞。

2.2 有限元模型

考慮到隧道縱向長度較長，計算斷面位于榆樹川隧道里程DK273+704處，該位置距離隧道進口30 m，隧道內空氣溫度可取為大氣溫度[15]。計算模型寬AD=BC=60 m，高AB=CD=80 m，隧道埋深15 m，計算模型如圖8所示，圍巖和襯砌的熱力學參數見表4。

圖8 隧道計算模型Fig.8 Numerical model of the tunnel

表4 圍巖和襯砌熱力學參數

2.3 邊界條件

溫度邊界：模型AB和CD為絕熱邊界；隧道襯砌邊界與隧道內空氣對流換熱，對流換熱系數為12.5 W/(m2·℃)。模型底部BC的熱流密度為0.01 W/m2。

水分邊界：考慮到隧道頂部的植被會減弱降雨和降雪對圍巖含水量的影響，同時在初期支護與二次襯砌之間施做防水層，因此，模型AB、BC、CD和AD邊界均為不透水邊界。

3 計算結果與分析

3.1 無保溫材料隧道溫度場分布

在隧道運營期間，隧道進口暴露在空氣之中，進口處圍巖的凍融狀態受空氣溫度的直接影響[20]。圖9為隧道運營20年后，襯砌支護背面最低溫度的分布。可以看出：隧道圍巖溫度沿隧道軸線呈對稱分布，最低溫度的最小值出現在隧道墻頂DL和DR處，最低溫度的最大值出現在仰拱下部G處。因此，將隧道的墻頂處DL和DR處作為隧道發生凍結的最不利位置。圖10為隧道運營20年內最不利位置處溫度隨時間變化。由圖可以看出，在隧道運營期，圍巖發生凍結的最不利位置處溫度隨時間呈周期性波動，并且最不利位置處的年最低溫度隨時間不斷降低，在隧道運營后第8年2月2日溫度降低至最低，并在隨后的年份內年最低溫度無明顯變化。因此，將隧道運營后第8年2月2日作為隧道內最不利位置發生凍結的最不利時間。

圖11和圖12分別給出了最不利凍結時間時，無保溫層隧道圍巖溫度和水分分布。由圖可以看出，在最不利凍結時間，圍巖0 ℃等溫線呈環向分布，并且在隧道圍巖處徑向形成厚度約為2.845 m的凍結圈。圍巖水分發生了明顯地重分布。沿著凍結圈，圍巖的未凍水含量和含冰量變化劇烈，在凍結圈內，靠近隧道襯砌支護的位置，圍巖含冰量急劇增大，并達到30%，遠遠超過了圍巖的初始含水量18%，這就是隧道發生凍脹破壞的原因。而在靠近0 ℃等溫線的位置處，圍巖未凍水含量為18%，含冰量為12%，圍巖的總含水量為30%，也超過了圍巖的初始含水量18%，這主要是由于在凍融循環作用下，圍巖內的水分向凍結鋒面發生了遷移。因此，對于嚴寒地區隧道需施做相應的保溫措施，使隧道襯砌不因圍巖發生凍結而破壞。

圖9 隧道運營第20年初期支護背面最低溫度分布Fig.9 The minimum temperature at the boundary of preliminary lining of the selected section in the 20th year

圖10 隧道運營20年內最不利位置溫度隨時間變化Fig.10 Temperature variation at the unfavorable position of tunnel during 20 years

圖11 最不利時間時無保溫層隧道圍巖溫度分布Fig.11 Temperature distribution of the selected section of tunnel without thermal insulation layer at the unfavorable time

圖12 最不利時間時無保溫隧道圍巖水分分布Fig.12 Moisture distribution of the selected section of tunnel without thermal insulation layer at the unfavorable time

3.2 大氣溫度對凍結圈的影響

大氣溫度處于負溫是隧道圍巖內水分發生凍結的原因，進而引發隧道襯砌的凍脹破壞。所以，對不同溫度條件下，凍融圈變化規律的研究必不可少，而溫度的變化主要體現在空氣年平均溫度的變化方面。圖13是空氣年平均溫度分別為1.82、2.32、3.32、4.32、5.32 ℃時隧道的溫度分布。可以看出，空氣溫度對隧道溫度場的分布影響巨大，隨著空氣年平均溫度的增大，隧道各位置處的凍結深度顯著降低，以隧道墻頂為例，空氣年平均溫度為1.82 ℃時，墻頂處的最大凍結深度為4.13 m；隨著空氣年平均溫度升高至2.32 ℃，墻頂處的最大凍結深度減小為3.575 m；當空氣年平均溫度升高至5.32 ℃，墻頂處的最大凍結深度減小為2.048 m。隧道最大凍結深度隨空氣年平均溫度的變化如圖14所示，對隧道最大凍結深度隨年平均溫度的變化擬合函數關系，得

圖13 不同空氣溫度條件下無保溫層隧道溫度分布Fig.13 Temperature distribution of the selected section of tunnel without thermal insulation layer under different atmospheric temperature

d=6.132 4Ta-0.651 4

(14)

式(14)中：d為隧道最大凍結深度，m；Ta為空氣年平均溫度，℃。

隧道內不同位置處的空氣溫度往往與隧道外空氣溫度不同，如果按隧道外空氣溫度確定隧道內不同深度處的最大凍結深度會存在誤差。運用擬合得到的隧道最大凍結深度與空氣年平均溫度的變化關系[式(14)]，只需知道隧道內不同深度處的年平均溫度，即可確定隧道內不同深度處的最大凍結深度，可以很好地可以解決上述存在的問題，為隧道內不同深度處最大凍結深度的確定提供依據。

圖14 隧道最大凍結深度隨空氣年平均溫度的變化Fig.14 Variation in the maximum frozen depth of tunnel under different atmospheric temperature

3.3 保溫材料對凍融圈的影響

在寒區隧道工程中，為了防止隧道發生凍脹破壞，常采用圍巖注漿，鋪設防水層，設置止水帶，設置保溫層等綜合措施，而作為嚴寒地區隧道防治凍害的有效措施，保溫層常常鋪設于初期支護與二次襯砌之間以達到防止圍巖發生凍脹的目的。現采用上述水-熱耦合模型，對鋪設導熱系數為0.025 W/(m·℃)，厚度為5 cm保溫材料的隧道運營8年后的溫度場進行分析。

圖15為最不利凍結時間時有保溫層隧道圍巖溫度分布。由圖可以看出，鋪設保溫層后，在最不利凍結時間，圍巖已無0 ℃等溫線存在，這說明保溫層的鋪設使嚴寒地區隧道圍巖凍結狀況得到了明顯地改善，降低了圍巖水分凍結的可能性，避免了襯砌凍害的產生。

圖15 最不利時間時有保溫層隧道圍巖溫度分布Fig.15 Temperature distribution of the selected section of tunnel with thermal insulation layer at the unfavorable time

以隧道墻頂位置處，取二次襯砌表面、初期支護表面和圍巖表面作為關鍵點，分析保溫層對隧道保溫效果的影響。圖16為隧道運營第8年，有、無保溫層隧道各關鍵點溫度隨時間的變化規律。由圖可以看出：隧道襯砌在無保溫層時，襯砌表面溫度受隧道內空氣溫度的直接影響，呈正、負溫度交替變化，在冷季形成的低溫凍結圈在暖季全部融化，如此寒暖交替，凍融反復。而鋪設保溫層后，雖隧道二次襯砌溫度還隨時間呈正、負溫度交替變化，但初期支護和圍巖表面溫度全年都為正溫，從而避免了圍巖水分的凍結。可見，保溫層的鋪設對嚴寒地區隧道防凍有重要的意義。

圖16 運營第8年無保溫層和有保溫層隧道關鍵點溫度變化Fig.16 Temperature variation at the classical position of tunnel without and with thermal insulation layer in the 8th year

影響保溫效果的一個重要因素是保溫材料的導熱系數。考慮到目前寒區隧道采用的保溫材料的導熱系數一般在0.03 W/(m·℃)以下[21]，下面選取五個典型的導熱系數λ作趨勢分析。圖17是隧道分別鋪設導熱系數為0.008、0.010、0.012、0.015、0.017 W/(m·℃)的保溫材料(鋪設厚度均為10 cm)時，隧道運行第7～10年時間內，初期支護與保溫層接觸位置溫度變化規律。可以看出，該位置溫度變化幅度隨保溫材料導熱系數λ的增加而增加，在隧道表面敷設厚度為10 cm，導熱系數為0.012 W/(m·℃)的保溫材料后，已經可以保證隧道不出現負溫的不利工況。

圖17 不同導熱系數條件下初期支護表面溫度變化規律Fig.17 Variation of temperature at the unfavorable position of preliminary lining in the condition of different thermal conductivity

表5為表面鋪設不同導熱系數的保溫材料，隧道運行第7～10年時間內各典型位置的最大凍深。由表可以看出，隧道墻頂的凍深均隨著保溫材料導熱系數的減小而減小，當采用厚度為10 cm，導熱系數為0.017 W/(m·℃)的保溫材料時，隧道墻頂最大凍深為53.45 cm；當采用導熱系數為0.015 W/(m·℃)的保溫材料后，墻頂處的凍深減小到了31.04 cm；減小了42%，而當采用導熱系數為0.012 W/(m·℃)的保溫材料時，墻頂已經不會出現負溫的情況，說明隧道防凍與保溫層導熱系數密切相關。

保溫材料對隧道的保溫效果不僅受到保溫材料類型的影響，還和保溫材料的厚度密切相關。圖18是隧道鋪設導熱系數為0.015 W/(m·℃)，厚度分別為8、10、13、14、15 cm的保溫材料時，隧道運行第7～10年時間內，初期支護與保溫層接觸位置溫度變化規律。可以看出，該位置處溫度變化振幅隨保溫層厚度的增加而減小，在隧道表面敷設導熱系數為0.015 W/(m·℃)厚度為13 cm的保溫材料已經能保證隧道不出現負溫的不利工況。

圖18 不同厚度條件下初期支護表面溫度變化規律Fig.18 Variation of temperature at the unfavorable position of preliminary lining in the condition of different thickness

表6為表面鋪設不同厚度的保溫材料，隧道運行第7～10年時間內各典型位置的最大凍深。由表可以看出，隧道墻頂的凍深隨著保溫材料厚度的增大而減小，當保溫層厚度為8 cm時，隧道墻頂最大凍深為46.7 cm；當保溫層厚度為10 cm時，凍深則減小到了31 cm，減小了33.6%；而當保溫層厚度為13 cm時，該位置則已經不會出現負溫的情況，說明隧道防凍與保溫層厚度密切相關。

表6 最大凍深與保溫層厚度的關系

4 結論

(1)本文建立的土體凍結過程中的水-熱耦合模型，考慮了水分遷移和冰-水相變潛熱對土體凍結的影響，并借助前人的實驗結果，通過數值計算的方法驗證耦合模型的準確性，計算結果與前人實驗結果十分吻合，驗證了本文建立的水-熱耦合模型在研究土體凍結過程中的有效性。

(2)隧道的墻頂處是隧道發生凍結的最不利位置，最不利位置處溫度隨時間呈周期性波動，并且最不利位置處的年最低溫度隨時間不斷降低，在隧道運營后第8年2月2日溫度降低至最低，以此作為隧道內最不利位置發生凍結的最不利時間。

(3)大氣溫度是隧道溫度場的分布決定因素，隨著大氣年平均溫度的增大，隧道各位置處的凍結深度顯著降低，隧道最大凍結深度與年平均溫度之間存在一定的關系。因此，可根據隧道內不同深度處的年平均溫度確定隧道內不同深度處的最大凍結深度。

(4)保溫層的鋪設使嚴寒地區隧道圍巖凍結狀況得到了明顯地改善，避免了襯砌凍害的產生。隧道的凍深隨保溫材料導熱系數的減小而減小，隨保溫層厚度的增加而減小。對于榆樹川隧道，敷設厚度為10 cm，導熱系數為0.012 W/(m·℃)的保溫材料，則墻頂已經不會出現0 ℃以下的情況；敷設導熱系數為0.015 W/(m·℃)，厚度為13 cm的保溫材料己經能保證該位置不出現負溫的情況。