王志杰,姜逸帆,李金宜,蔣新政,周 平, *,雷飛亞
(1. 西南交通大學 交通隧道工程教育部重點實驗室,四川 成都 610031;2. 西南交通大學土木工程學院,四川 成都 610031)
隨著我國基礎(chǔ)設(shè)施建設(shè)的不斷發(fā)展,大量公路隧道在高海拔、高寒等惡劣條件下進行修建。寒區(qū)隧道通常受到季節(jié)性凍融、凍脹作用的影響,給隧道的安全運營帶來了極大的風險。目前,為了解決隧道凍害問題,國內(nèi)較為常用的抗防凍措施為鋪設(shè)保溫層。相關(guān)研究也證實了保溫層的有效性[1]。
受空氣流動影響,寒冷空氣對隧道的影響會延伸至洞內(nèi)一定深度。為保證隧道在運營期間不出現(xiàn)凍害,需要對寒區(qū)隧道的縱向溫度場分布和保溫層的鋪設(shè)長度展開研究。
現(xiàn)場實測作為最直接的研究方法,被廣泛應(yīng)用于溫度場研究。Wheeler等[2]通過現(xiàn)場測試研究了隧道風速與溫度場的均勻性,并得出了空氣速度對隧道溫度場的影響規(guī)律。Yan等[3]通過對德格隧道進行監(jiān)測,結(jié)合數(shù)值模擬等研究方法,得出隧道內(nèi)溫度場的時空演化規(guī)律。高焱等[4]通過現(xiàn)場實測,明確了綏陽隧道的縱向溫度場及風速分布狀況。
為了更深入地研究溫度場分布規(guī)律,國內(nèi)外學者展開了大量的理論研究。Kawamura等[5]根據(jù)疊加原理和能量守恒定律,建立了寒區(qū)隧道洞內(nèi)空氣流體和圍巖固體對流換熱的計算模型,并給出該計算模型的解析解。Krarti等[6]通過對某地下風洞進行研究,得到了地下洞室的溫度場解析解,并用溫度的年振幅和年平均氣溫2個參數(shù)對溫度場特性進行了表征。何春雄等[7]建立了隧道內(nèi)空氣與圍巖對流換熱及固體導(dǎo)熱的綜合模型,分析預(yù)報了祁連山區(qū)大坂山隧道開通運營后的圍巖溫度分布情況。韓躍杰等[8]通過建立隧道縱向洞內(nèi)空氣與洞壁的氣固耦合傳熱模型,計算得到了隧道的縱向溫度場分布情況。
隨著計算機軟件的普及,通過軟件模擬隧道的溫度場也成為了更多學者的選擇。鄭陽等[9]通過Fluent與Ansys的假耦合計算,建立隧道圍巖、襯砌及洞內(nèi)氣體的流固對流換熱模型,得出了隧道圍巖的溫度場分布情況,并驗證了高寒隧道保溫層的必要性。賈輝[10]通過CFD數(shù)值模擬計算,分析了氣象因素對寒區(qū)隧道溫度場的影響規(guī)律及相關(guān)性。
對于如何合適地選取保溫層的鋪設(shè)長度,部分學者展開了相應(yīng)研究。孫文昊[11]根據(jù)現(xiàn)場測溫結(jié)果,得到洞內(nèi)空氣溫度的分布特征,綜合考慮保溫層的隔熱效應(yīng)修正和邊界效應(yīng)修正,得到隧道保溫層鋪設(shè)長度的計算公式。夏才初等[12]基于由隧道進、出口氣象條件及隧道地形條件求解的溫度場解析解,提出了一種保溫層鋪設(shè)長度計算方法。高焱[13]統(tǒng)計分析了大量寒區(qū)隧道的凍害情況,對設(shè)防長度與洞口氣溫關(guān)系進行擬合分析,得到了設(shè)防長度的統(tǒng)計規(guī)律。
目前已經(jīng)有不少學者展開了寒區(qū)隧道縱向溫度場的研究,且結(jié)果較為成熟,但是如何計算保溫層設(shè)防長度,還處于一個不斷探索的階段。隨著國家公路建設(shè)的發(fā)展,由于地形和選線的原因,不可避免地出現(xiàn)大曲率曲線隧道,甚至螺旋隧道,現(xiàn)有的研究沒有涉及曲線隧道的溫度場計算及對應(yīng)的保溫設(shè)計需求。為了探究曲線隧道的保溫層鋪設(shè)長度計算方法,本文采用文獻調(diào)研、理論推導(dǎo)、數(shù)值模擬的方法,對曲線螺旋隧道的溫度場展開研究,結(jié)合理論推導(dǎo)提出一種較為可靠的曲線寒區(qū)隧道保溫層鋪設(shè)長度計算方法。
金家莊特長螺旋隧道為河北省政府和交通運輸部“一號工程”——延崇高速公路的重點控制性工程,其入口位于赤城縣炮梁鄉(xiāng)金家莊村西北方向,出口位于崇禮區(qū)棋盤梁村,左幅長4 228 m,右幅長4 104 m。隧道設(shè)斜井1處,長750.95 m,按雙向4車道高速公路標準建設(shè),最大行車速度為100 km/h。建設(shè)螺旋隧道,能解決隧道進出口的大高差和路線縱坡太大的問題,保證車輛安全行駛。隧道設(shè)計曲率半徑達到890 m,隧道平面示意圖見圖1。
圖1 金家莊特長公路隧道平面示意圖
隧址區(qū)位于大陸性季風氣候中溫帶亞干旱區(qū),四季分明,冬季寒冷漫長,降雪量較少。全年日照時間長、溫差大,風向以西北風和靜風為主。年均降水量在424 mm 左右,且分布不均,全年無霜期平均為115.9 d。年平均氣溫為5.5 ℃,年最髙氣溫為39.4 ℃,最低氣溫為-28.2 ℃,最大積雪深度為9 cm,最大凍土深度為162 cm,無霜期為145 d。
隧道是一個復(fù)雜的結(jié)構(gòu)體,隧道內(nèi)的空氣流動形式也同樣復(fù)雜。為了得到隧道保溫層鋪設(shè)長度的計算公式,進行如下假設(shè):
1)研究對象達到平衡狀態(tài),即為穩(wěn)態(tài)問題;
2)只考慮徑向傳熱,忽略隧道空氣縱向?qū)α鲹Q熱;
3)所有材料均為各向同性,且導(dǎo)熱系數(shù)、比熱容、密度不隨溫度變化;
4)不考慮二次襯砌與初襯、襯砌與圍巖之間的接觸熱阻,接觸邊界處滿足溫度和熱流量相等的連續(xù)條件;
5)將復(fù)雜的傳熱問題簡化為帶有對流換熱邊界的三維多層圓筒壁的熱傳導(dǎo)問題,將隧道、襯砌及圍巖斷面簡化為圓形或圓環(huán)。
多層圓筒傳熱模型作為求解傳熱模型的經(jīng)典模型,已經(jīng)被引入隔熱層厚度計算[14]。對于厚度計算求解,使用一維傳熱模型即可; 而對于鋪設(shè)長度的計算,需要引入三維模型,考慮縱向方向的溫度傳遞。因此,本文采用三維的多層圓筒傳熱模型。
簡化后的傳熱模型如圖2所示。
r1、r2、r3、r4分別為二次襯砌與空氣交界面、初期支護與二次襯砌交界面、初期支護與圍巖交界面、圍巖恒溫邊界距圓筒中心的距離; λ1為二次襯砌的導(dǎo)熱系數(shù); λ2為初期支護的導(dǎo)熱系數(shù); λ3為圍巖的導(dǎo)熱系數(shù); t1(z)為進洞深度z(m)的空氣溫度; t2(z)為進洞深度z(m)的襯砌背后圍巖溫度; tw為進洞深度z(m)的深處圍巖恒溫溫度。
1)根據(jù)已有的研究,多層圓筒傳熱模型徑向的傳熱基本方程為:
(1)
式中:Q為熱通量;R為圓筒內(nèi)外兩側(cè)的熱阻; Δt為圓筒內(nèi)外兩側(cè)的溫差。
(2)
式中:Rh為單位長度空氣圍巖換熱等效熱阻;R1為單位長度二次襯砌支護等效熱阻;R2為單位長度初期支護等效熱阻;R3為單位長度圍巖等效熱阻。
(3)
式中hc為對流換熱系數(shù)。
(4)
(5)
(6)
2)為了計算縱向的溫度傳遞,本文做如下處理。如圖3所示,取空氣微元dz,空氣吸收圍巖帶來的能量為Φ,空氣溫度由t1(z)變?yōu)閠1(z+Δz),空氣微元dz溫度變化的控制方程為:
(7)
式中:ρ為空氣密度;c為空氣比熱容;A為隧道斷面面積。
圖3 空氣微元吸熱示意圖
根據(jù)圓筒傳熱方程,圍巖傳向空氣的熱量即為空氣吸收的熱量,所以微元dz吸收的熱量
(8)
式中v為空氣速度。
3)邊界條件為:
(9)
式中t0為洞口空氣溫度。
4)保溫層鋪設(shè)長度計算:
聯(lián)立式(7)—(9)可得
(10)
求解微分方程得
(11)
則二次襯砌背后溫度
(12)
由式(12)可以得到保溫層鋪設(shè)長度
(13)
表面對流換熱系數(shù)的數(shù)值與換熱過程中流體的物理性質(zhì)、換熱表面的形狀、部位以及流體的流速等都有密切關(guān)系。物體表面附近的流體的流速越大,其表面對流換熱系數(shù)也越大。計算傳熱系數(shù)的方法主要有試驗求解法、數(shù)學分析解法和數(shù)值分析解法。
國內(nèi)已有學者對混凝土的對流換熱系數(shù)進行了試驗求解。張建榮等[15]利用風洞探究了不同風速下的對流換熱系數(shù)hc,總結(jié)出如下關(guān)系式:
hc=3.06v+4.11 。
(14)
不同曲率的隧道在環(huán)境溫度以及初始地溫、斷面大小等相同的情況下,其對流換熱的差異可以簡化認為是由空氣速度的差異造成。
利用CFD軟件Fluent探究不同曲率隧道的空氣速度影響。
3.2.1 模型的建立
采用控制變量法,保持隧道總長度(4 220 m)不變,為了使隧道圍巖外邊界溫度保持為原始巖溫,所建模型的橫斷面尺寸應(yīng)超出調(diào)熱圈范圍[16],計算模型橫斷面采用邊長40 m的正方形區(qū)域,隧道凈空面積為85.42 m2。本文模型的初始溫度場通過UDF賦予初始地溫實現(xiàn),所賦予的初始地溫考慮了洞外空氣溫度以及地熱的影響,因此,賦予初始地溫后的模型不再考慮邊界溫度的影響,采用絕熱邊界。通過改變隧道曲率半徑,設(shè)置曲率半徑r為900、1 800、2 700、3 600、4 500、5 400、6 300 m的螺旋隧道及同長度直線隧道共8個工況。隧道模型如圖4所示。
(a) 工況1: r=900 m螺旋隧道(b) 工況2: r=1 800 m螺旋隧道
(c) 工況3: r=2 700 m螺旋隧道(d) 工況4: r=3 600 m螺旋隧道
(e) 工況5: r=4 500 m螺旋隧道(f) 工況6: r=5 400 m螺旋隧道
(g) 工況7: r=6 300 m螺旋隧道(h) 工況8: 直線隧道
(i) Fluent計算模型 (單位: m)
3.2.2 初始地溫的計算
由于隧道從進口到出口各處的埋深不同,造成了隧道全線不同位置的初始地溫存在差異,因此,準確地模擬出初始地溫可以提高計算的準確性。本文采用的初始地溫處理流程如圖5所示。
圖5 初始地溫處理流程
根據(jù)氣象部門提供的隧址區(qū)大氣溫度數(shù)據(jù)及現(xiàn)場氣象站監(jiān)測數(shù)據(jù),繪制日平均溫度變化時程曲線。通過三角函數(shù)對隧址區(qū)日平均溫度時程變化曲線進行擬合,擬合函數(shù)如下:
Temp(t′)=5.5-20.5sin(2πt′/365) 。
(15)
式中: Temp(t′)為隧址區(qū)洞口處大氣溫度擬合函數(shù),℃;t′為時間,d。
確定初始溫度場時,模型的初始溫度考慮地熱,將模型在年溫度周期循環(huán)200年后的溫度場作為金家莊特長公路隧道初始地溫場,根據(jù)地質(zhì)調(diào)查圖紙,利用Ansys建立平面二維模型。模型的邊界條件加載如下: 模型上邊界與大氣接觸,邊界條件設(shè)定為隨時間動態(tài)變化的溫度變化荷載Temp(t′); 下邊界與深處的地層存在換熱,設(shè)為熱流密度邊界,根據(jù)相關(guān)文獻熱流密度取q=0.06 W/m2; 模型左、右邊界由于與接觸的圍巖處于同一埋深,受外界環(huán)境的影響相同,因此,假設(shè)沒有熱量的傳遞,設(shè)為絕熱邊界。金家莊特長公路隧道初始地溫場計算模型及邊界條件如圖6所示。
(a) 數(shù)值計算模型
(b) 數(shù)值計算邊界
圖7為計算得到的200年后金家莊特長公路隧道冬季(1月1日)溫度場分布云圖。
圖7 金家莊特長公路隧道冬季溫度場分布云圖(單位: ℃)
圖8示出金家莊特長公路隧道中間深埋段冬季(1月1日)圍巖溫度場分布曲線。由圖可以看出,在隧道的中間深埋段,圍巖溫度隨里程變化較小(12.5~13.5 ℃),圍巖溫度分布曲線表現(xiàn)出與地形相似的幾何特征,即中間高、兩端低的分布趨勢。
圖8 金家莊特長公路隧道中間深埋段冬季圍巖溫度場分布曲線
圖9示出金家莊特長公路隧道進出口淺埋段冬季(1月1日)圍巖溫度場分布曲線。由圖可以看出: 隧道洞口段(100 m范圍)溫度分布受外界氣象溫度影響明顯,隨著進洞深度的增加,隧道襯砌位置的圍巖溫度出現(xiàn)迅速增大的趨勢,變化幅度接近6 ℃,在建模計算過程中不可忽略。為了得到更加符合實際的圍巖溫度初始條件,通過函數(shù)擬合的方式得到隧道進出口段的圍巖溫度隨進洞深度x變化函數(shù),擬合公式為:
(16)
將擬合的溫度場采用UDF編程,賦予Fluent簡化模型。隧道入口和出口段的初始溫度場如圖10所示,基本符合實際初始地溫。
(a) 隧道進口圍巖初始溫度
(b) 隧道出口圍巖初始溫度
(a) 隧道入口 (b) 隧道出口
3.2.3 計算參數(shù)
根據(jù)隧址區(qū)日平均溫度時程變化曲線的擬合函數(shù),設(shè)置洞口氣溫為隧址區(qū)最冷月平均氣溫-15 ℃; 根據(jù)隧道自然風要求,進口風壓設(shè)為28.49 Pa(對應(yīng)風速2 m/s),進口風速取單向持續(xù)通風60 d計算。熱力學計算參數(shù)如表1所示。
表1 熱力學計算參數(shù)
3.2.4 計算結(jié)果分析
不同工況(隧道曲率半徑不同)的空氣速度斷面圖見圖11。
通過對比不同曲率半徑的空氣流速特點,可以得出以下結(jié)論:
1)由于曲率的影響,空氣在隧道內(nèi)流動中會受到離心力的作用,向曲率隧道的外側(cè)偏移,擠壓氣體,使外側(cè)的氣體擁有更高的流動速度,而內(nèi)側(cè)的氣體流速相對減小。
2)曲率與速度的偏移程度正相關(guān),曲率越大速度的偏移現(xiàn)象越明顯。
根據(jù)不同曲率半徑螺旋隧道及直線隧道的計算結(jié)果,在進洞50、150、300、500、800、1 100、1 500、2 000 m深度位置共設(shè)置8個監(jiān)測斷面(編號1—8),布置從隧道左側(cè)到右側(cè)的測線,可以得到不同曲率半徑隧道在不同進洞深度的斷面空氣速度分布曲線,如圖12所示。
根據(jù)圖12可以得到以下結(jié)論:
1)速度的偏移現(xiàn)象需要在空氣流進入隧道一定深度后才會出現(xiàn)。一般是在第2個斷面即洞深150 m后才出現(xiàn)偏移,且偏移發(fā)生后會迅速達到穩(wěn)定狀態(tài),即隧道中的空氣流動將保持持續(xù)的偏移狀態(tài)。
2)由于混凝土壁面的摩擦影響,空氣在邊界層位置的速度明顯小于非邊界位置的速度。偏移作用的擠壓加速效應(yīng)在邊界層上依然適用,依然呈現(xiàn)曲率越大、速度偏移越大的規(guī)律。
由于與混凝土發(fā)生對流換熱的空氣處于邊界層,因此,研究邊界層的空氣速度更有意義。其空氣速度越大,設(shè)防長度越長。圖13示出不同曲率半徑隧道的邊界層空氣速度曲線。
(a) 工況1: r=900 m螺旋隧道 (b) 工況2: r=1 800 m螺旋隧道
(c) 工況3: r=2 700 m螺旋隧道 (d) 工況4: r=3 600 m螺旋隧道
(e) 工況5: r=4 500 m螺旋隧道 (f) 工況6: r=5 400 m螺旋隧道
(g) 工況7: r=6 300 m螺旋隧道 (h) 工況8: 直線隧道
(a) r=900 m螺旋隧道 (b) r=1 800 m螺旋隧道 (c) r=2 700 m螺旋隧道
(d) r=3 600 m螺旋隧道 (e) r=4 500 m螺旋隧道 (f) r=5 400 m螺旋隧道
(g) r=6 300 m螺旋隧道 (h) 直線隧道
圖13 不同曲率半徑隧道的邊界層空氣速度曲線
對圖13曲線進行擬合,擬合度為0.996,擬合效果良好,擬合公式為:
(17)
式中r為隧道曲率半徑。
將得到的隧道曲率半徑與邊界層空氣速度的擬合公式(17)以及混凝土對流換熱系數(shù)公式(14)代入保溫層鋪設(shè)長度計算公式(13)中,得到:
(18)
根據(jù)曲線隧道保溫層鋪設(shè)長度計算公式(18)可以得到不同曲率半徑隧道的保溫層鋪設(shè)長度,如表2所示。
表2 不同曲率半徑隧道的保溫層鋪設(shè)長度
不同曲率半徑的隧道二次襯砌背后溫度數(shù)值模擬結(jié)果如圖14所示,二次襯砌背后凍深長度模擬計算結(jié)果如表3所示。
圖14 不同曲率半徑的隧道二次襯砌背后溫度曲線數(shù)值模擬結(jié)果
表3 二次襯砌背后凍深長度模擬計算結(jié)果
將數(shù)值模擬計算結(jié)果與公式計算結(jié)果進行對比,如圖15所示。
圖15 保溫層鋪設(shè)長度數(shù)值模擬計算與公式計算結(jié)果對比
可以看出,曲線隧道的保溫層鋪設(shè)長度公式計算結(jié)果與數(shù)值模擬計算結(jié)果對比,呈現(xiàn)相同的變化規(guī)律,且鋪設(shè)長度隨曲率半徑的變化規(guī)律大致相同,證明了計算公式的可靠性。
本文依托金家莊隧道,通過圓筒傳熱模型推導(dǎo)保溫層鋪設(shè)的理論計算公式,并利用CFD軟件探討了不同曲率半徑隧道溫度場分布規(guī)律以及保溫層長度鋪設(shè)規(guī)律,得出了曲率影響隧道溫度場的關(guān)鍵因素,結(jié)合理論推導(dǎo)得出了曲線隧道保溫層鋪設(shè)長度計算公式。
1)根據(jù)圓筒傳熱模型推導(dǎo)出的保溫層鋪設(shè)長度計算公式可以發(fā)現(xiàn),不同的隧道曲率對鋪設(shè)長度的影響主要反映在對邊界層空氣速度的影響。不同的空氣速度影響了空氣與混凝土襯砌的對流換熱系數(shù)以及空氣受熱加熱時間,進而影響了溫度的傳遞。
2)通過數(shù)值模擬,發(fā)現(xiàn)曲線寒區(qū)隧道受曲率的影響,隧道內(nèi)斷面空氣峰值速度在離心力的作用下向隧道外側(cè)偏移,偏移程度與曲率半徑正相關(guān)。空氣的偏移對外側(cè)的邊界層產(chǎn)生擠壓效應(yīng),加速了外側(cè)邊界層的空氣流動速度。偏移效應(yīng)一般在隧道進洞深度150 m出現(xiàn)并迅速達到穩(wěn)定。
3)通過溫度場分析,曲線隧道的溫度場分布沿隧道軸線呈非對稱分布,外側(cè)邊界層空氣在離心力作用下加快了對流換熱,因此,外側(cè)的凍深長度大于內(nèi)側(cè),在實際設(shè)防時應(yīng)重點關(guān)注外側(cè)的襯砌背后溫度變化情況,且隧道曲率越大該現(xiàn)象越明顯。
4)根據(jù)圓筒傳熱模型推導(dǎo)出的保溫層鋪設(shè)長度計算公式,是一種簡化計算公式。該計算公式可用于計算不同曲率半徑隧道的保溫層鋪設(shè)長度,與數(shù)值模擬結(jié)果對比,呈現(xiàn)相同的變化規(guī)律,但計算結(jié)果較為保守,能夠滿足對抗防凍的鋪設(shè)長度要求,在實際應(yīng)用時可考慮20%的折減。后續(xù)建議耦合多種影響因素對保溫層的鋪設(shè)長度進行研究。