基于圓形斷面的隧道溫度場有限差分計算模型

周小涵， 曾艷華， 白 赟， 魏英杰, 何 省

(1． 西南交通大學交通隧道工程教育部重點實驗室， 四川 成都 610031； 2． 中鐵大橋勘測設計院集團有限公司， 湖北 武漢 430050； 3． 成都建筑工程集團總公司， 四川 成都 610014；4． 中國中鐵二院工程集團有限責任公司， 四川 成都 610031)

?

基于圓形斷面的隧道溫度場有限差分計算模型

周小涵1，2， 曾艷華1,*， 白 赟1， 魏英杰3, 何 省4

(1． 西南交通大學交通隧道工程教育部重點實驗室， 四川 成都 610031； 2． 中鐵大橋勘測設計院集團有限公司， 湖北 武漢 430050； 3． 成都建筑工程集團總公司， 四川 成都 610014；4． 中國中鐵二院工程集團有限責任公司， 四川 成都 610031)

除了少數圓形隧道外，大部分山嶺交通隧道斷面采用馬蹄形或端墻式等形狀。在隧道溫度場的預測中，圓形斷面模型能否代替馬蹄形等實際隧道模型，其適應性值得研究。應用基于空氣-襯砌-圍巖的對流-導熱耦合作用控制方程的有限差分方法，建立圓形斷面模型對東北寒區馬蹄形隧道溫度場進行計算，并與現場實測溫度場進行對比。結果表明： 1) 圓形斷面隧道模型有限差分計算方法克服了通用有限元軟件建模復雜、對硬件要求高的弊端，考慮了隧道內風流速度和入口風流溫度的影響，在隧道溫度場的預測計算中能夠滿足工程使用要求。2)隧道內風流速度和入口風流溫度對隧道溫度場影響較大。本文算例中，入口風流溫度每升高10 ℃，二次襯砌表面溫度升高約7.2 ℃，增幅均勻；從1～5 m/s,洞內風流速度每增大1 m/s,二次襯砌表面溫度降低的幅度為6.6、2.7、1.5、0.9 ℃，降幅越來越小。

馬蹄形隧道； 圓形隧道； 溫度場； 有限差分； 現場實測； 風流溫度； 風流速度

0 引言

準確預測隧道溫度場在隧道貫通后的變化規律是寒區和高地溫隧道設計的前提，國內外學者對隧道溫度場展開過多方面的研究。

在隧道結構的熱傳導方面： R.N.Pande等[1]推導了兩相介質的整合熱傳導率取值方法；J.A.Krumhansl[2]討論了有序和無序固體介質的熱傳導問題；張耀等[3]考慮正弦曲線規律變化的對流換熱邊界條件，建立了復合結構的圓形隧道熱傳導方程。

在理論求解方面： 賴遠明等[4]應用無量綱量和攝動技術，給出了圓形隧道凍結過程的近似解析解；S.Singh等[5-6]利用疊加原理和分離變量方法得到了圓形斷面隧道在考慮對流換熱邊界條件下的瞬態溫度場解析解；張國柱等[7]利用疊加原理及貝塞爾特征函數的正交及展開定理,得到了寒區隧道圍巖徑向溫度的理論解;J.L.Battaglia等[8]進行了傳熱模型的拉普拉斯變換反演研究；邵珠山等[9]利用無量綱化和微分方程級數求解法,得到了溫度場、位移場及應力場的熱彈性理論解；鄭陽等[10]通過Fluent與Ansys的假耦合計算,建立了馬蹄形隧道圍巖、襯砌及洞內氣體的流固對流換熱模型。

在保溫隔熱層計算方法方面： K.Okada[11]提出了寒區隧道隔熱層的一種設計方法；陳建勛[12]推導了圓形隧道模型防凍隔溫層厚度的計算方法；夏才初等[13]采用理論分析和數值模擬計算等手段，研究了圓形隧道不同類型凍土段隔熱(保溫)層鋪設厚度計算方法。

在溫度實測方面： 陳建勛等[14]對隧道拱頂、拱腰邊墻和路面4個部位的溫度進行1.5年的測試，總結了其溫度變化規律；郝飛等[15]利用Ansys軟件對馬蹄形寒區隧道進行了模擬分析,將計算的理論值與實測值進行了對比。

除了少數圓形隧道外，大部分公路和鐵路隧道斷面采用的是馬蹄形或端墻式等形狀，圓形斷面模型能否代替馬蹄形等實際隧道模型，其適應性值得研究。在當前的實際應用中，若使用通用計算軟件(Fluent等)，可采用隧道實際斷面模型；但在實際建模中，由于該類軟件無法同時考慮對流換熱、圍巖導熱和空氣傳熱等問題，隧道圍巖原始溫度等邊界條件設置困難，模型網格量大，對電腦硬件要求高，工程適應性較差。

文獻[16]建立了基于空氣-襯砌-圍巖的對流-導熱耦合控制方程的有限差分計算模型，本文旨在通過研究以圓形斷面代替馬蹄形斷面計算隧道溫度場的適應性，以驗證基于空氣-襯砌-圍巖的對流-導熱耦合控制方程的有限差分計算模型在隧道溫度場預測中的可行性，同時重點分析隧道內風流速度和隧道入口風流溫度2大因素對隧道溫度場的影響。

1 有限差分控制方程

文獻[16-17]在已有研究中基于空氣-襯砌-圍巖的對流-導熱耦合控制方程，用有限差分法實現了計算，同時考慮了隧道入口風流溫度和洞內風流速度對隧道溫度場的影響。本文列出圍巖和氣流溫度場主要控制方程，其他方程不再贅述。

1.1 圓形有限差分模型控制方程

在流體力學中，把水力半徑相等的圓管直徑定義為非圓管的當量直徑

R=A/L=0.25d。

(1)

式中:A=(π×d2)/4，為截面積;L=π×d,為圓周長度。與圓形斷面相比，非圓形斷面的當量直徑也可以用4倍的水力半徑表示。

隧道縱向瞬態溫度場圍巖的控制方程(柱坐標系下)可以表示為

(2)

邊界條件

將隧道內氣流視為無黏性、不可壓縮的穩定流體，根據隧道內縱向氣流的質量守恒、動量守恒和能量守恒定律，得到隧道縱向瞬態溫度場氣流部分的控制方程

(5)

邊界條件為

tf(x,τ)=tb(τ=0)。

(6)

式(2)—(6)中: t是襯砌或圍巖內部任一點的溫度,℃； r是襯砌或圍巖中任一點到隧道中心的距離，m； τ是時間參數,s； ρ是空氣的密度，kg/m3； a為圍巖的導溫系數； cp是空氣的比熱容,J/(kg·℃)； tf是風流溫度,℃； tb是隧道壁面溫度,℃； tw是圍巖溫度,℃； h是對流換熱系數,W/(m2·℃)； U是隧道斷面周長,m； qs是隧道內熱源的發熱量,W； v是風流速度,m/s。

1.2 有限差分相關參數

為了在實際隧道溫度場計算中應用，必須根據隧道的實際結構對其空間域進行離散處理。

沿徑向，將二次襯砌、初期支護分成5等分，并將圍巖劃分為若干個子區域，每個子區域間距為1m，節點1為對流邊界，節點6和節點11為復合材料傳熱邊界，其余點均為內部點。計算時沿縱向每30m設置一個節點。橫斷面內差分節點劃分如圖1所示。縱斷面內差分節點劃分如圖2所示。

圖1 橫斷面內差分節點劃分

Fig. 1 Partition of primary support, secondary lining and surrounding rock

圖2 縱斷面內差分節點劃分

用偏微分方程的替代法建立隧道凈空與襯砌之間、襯砌內部、襯砌與圍巖之間、圍巖內部的傳熱及隧道氣流傳熱的中心有限差分方程，結合能量守恒定律建立三維瞬態傳熱的異步長顯示差分格式。利用差分格式穩定性判據式，確定時間步長Δτ與隧道縱向空間步長Δx，編制程序進行計算。

2 寒區隧道溫度實測

2015年1月，對東北牡綏線鐵路改造工程在建長度為6 170 m的綏陽隧道(于2013年9月貫通)進行溫度實測。根據隧道長度在綏陽隧道設置了13個測試斷面，如圖3所示。隧道二次襯砌表面溫度測試點由拱頂到軌道面依次布設，如圖4所示。使用紅外線測溫儀GM1150對隧道斷面二次襯砌表面的4個測點進行溫度實測和記錄分析。使用機械風表對隧道內風流速度和風向進行了測試，測試結果表明，綏陽隧道內的自然風流速度為3.8～4.0 m/s，風向都是自小里程端吹向大里程端。

3 溫度場計算值與實測值對比

綏陽隧道橫斷面凈空面積約為94 m2，橫斷面周長約為37.5 m，取其等效水力半徑為2.5 m。對綏陽隧道所在地區從隧道貫通到溫度實測期間的溫度進行統計分析，得出隧道入口風流溫度tf(℃)隨時間i(d)的關系

tf=2.5+18.5×sin (2×π×i /365)。

(7)

圖3 綏陽隧道溫度測試斷面隨縱向分布示意圖

圖4 隧道二次襯砌表面溫度測試點橫斷面分布圖

使用tf作為隧道入口風流溫度邊界條件進行計算，計算風流速度取隧道內實測平均風流速度3.8 m/s，方向由小里程端吹向大里程端。隧道原始巖溫按照地勘報告中埋深和地溫梯度3.0 ℃/100 m選取。根據風流速度與對流換熱系數的關系[18]，本次計算對流換熱系數取為21 W/(m2·℃)，其他介質材料計算參數取值見表1。由于缺少實驗參數，本文計算不考慮相變問題。

應用式(2)—(6)的有限差分控制方程，編程計算得出綏陽隧道由2013年9月貫通至2015年1月隧道襯砌表面溫度場分布情況。將綏陽隧道2015年1月7—13日各襯砌表面測點的實測溫度平均值和數值計算結果中相應時段的襯砌表面平均溫度進行比較分析，如圖5所示。

表1 介質材料參數

圖5 綏陽隧道襯砌表面溫度實測值和圓形斷面模型計算值

Fig. 5 Comparison between measured lining surface temperature of Suiyang Tunnel and circular cross-section model calculated results

通過現場溫度測試數據可以看出，隧道拱頂二次襯砌表面溫度最高而軌道面測點溫度最低。從拱頂到軌道面，4個測點溫度逐漸降低。各測點溫度平均值與測點2和測點3的溫度接近，測點平均溫度比拱頂測點溫度低約0.7 ℃，比軌道面測點溫度高約0.4 ℃。

對比有限差分計算結果和實測溫度結果可知，兩者隨隧道縱向分布趨勢一致且差值較小。隧道出入口段計算溫度和實測溫度有較大差異，主要是因為洞口處圍巖受外界大氣影響顯著，而該因素在模型計算中考慮不足以及計算中沒有考慮相變所致。在洞身段，數值計算襯砌平均溫度高于實際測點平均溫度約1.9 ℃，滿足工程使用要求。

4 風流溫度和風流速度對溫度場的影響

基于本文有限差分算法，利用綏陽隧道實際模型，重點研究隧道內風流溫度和風流速度對溫度場的影響。為了便于研究，分2種情況： 1)保持隧道內風流速度3.8 m/s 不變，隧道入口風流溫度采用-20、-10、0、10、20 ℃ 5個恒定值；2)保持入口風流溫度-20 ℃不變，隧道內風流速度分別采用1、2、3、4、5 m/s 5個恒定值。選擇隧道橫斷面DK493+529的計算結果進行分析。不同風流溫度和風流速度條件下隧道二次襯砌表面溫度如圖6和圖7所示。

圖6 不同風流溫度時綏陽隧道襯砌表面溫度(DK493+529)

Fig. 6 Lining surface temperature of Suiyang Tunnel under different wind flow temperatures at cross-section DK493+529

圖7 不同洞內風流速度條件下綏陽隧道襯砌表面溫度(DK493+529)

Fig. 7 Lining surface temperature of Suiyang Tunnel under different wind flow velocities at cross-section DK493+529

從圖6可知，隧道入口風流溫度對隧道襯砌溫度影響較大。本文算例中，計算時間為30 d時，風流溫度每升高10 ℃，二次襯砌表面溫度升高約7.2 ℃，增幅均勻。

從圖7可知，隧道洞內風流速度對隧道襯砌溫度影響較大。當有低溫風流對隧道進行降溫作用時，隨著風流速度的增大，二次襯砌表面溫度越來越低，但是其降低的幅度越來越小。計算時間為30 d時，從1～5 m/s,洞內風流速度每增大1 m/s,二次襯砌表面溫度降低的幅度為6.6、2.7、1.5、0.9 ℃。

5 結論與討論

應用基于空氣-襯砌-圍巖的對流-導熱耦合作用控制方程的有限差分方法，建立圓形斷面模型對東北寒區馬蹄形隧道溫度場進行計算，并與實測溫度場進行對比，得出如下結論。

1)根據現場溫度測試數據可知，隧道拱頂二次襯砌表面的溫度最高而軌道面測點溫度最低。各測點溫度平均值與拱腰的溫度接近。對比寒區隧道溫度場有限差分計算結果和實測結果可知，兩者隨隧道縱向分布趨勢一致且差值較小。圓形斷面模型有限差分計算方法克服了通用有限元軟件無法同時考慮對流換熱和傳熱問題、建模復雜和對硬件要求高等弊端，在隧道溫度場長期分布規律的預測計算中易于使用。

2)隧道入口風流溫度和洞內風流速度對隧道溫度場有較大影響。本文算例中，入口風流溫度每升高10 ℃，二次襯砌表面溫度升高約7.2 ℃，增幅均勻； 從1～5 m/s,洞內風流速度每增大1 m/s,二次襯砌表面溫度降低的幅度為6.6、2.7、1.5、0.9 ℃，降幅越來越小。

3)在收集隧道當地氣候資料而得到預測隧道入口風流溫度的條件下，本文計算模型可以用來計算不同隧道內風流速度條件下的溫度場分布，為具體設計提供幫助。本方法對凍土(季節性凍土和永久性凍土)隧道的溫度場及隔熱層厚度預測的適應性值得進一步研究。

[1] Pande R N, Kumar V, Chaudhary D R. Thermal conduction in a homogeneous two-phase system[J]. Pramana, 1984, 22(1): 63-70.

[2] Krumhansl J A. On the general physics of heat conduction in ordered and disordered solids[J]. Springer Berlin Heidelberg, 1993, 112: 3-8.

[3] 張耀，何樹生，李靖波. 寒區有隔熱層的圓形隧道溫度場解析解[J]. 冰川凍土, 2009, 31(1): 113-118. (ZHANG Yao，HE Shusheng，LI Jingbo. Analytic solutions for the temperature fields of a circular tunnel with insulation layer in cold region[J]. Journal of Glaciology and Geocryology，2009，31(1)： 113-118. (in Chinese))

[4] 賴遠明，喻文兵，吳紫汪，等. 寒區圓形截面隧道溫度場的解析解[J]. 冰川凍土, 2001, 23(2): 126-130. (LAI Yuanming, YU Wenbing, WU Ziwang, et al. Approximate analytical solution for the temperature fields of a circular tunnel in cold regions[J]. Journal of Glaciology and Geocryology, 2001, 23(2): 126-130. (in Chinese))

[5] Singh S，Jain P K，Rizwan-uddin. Analytical solution to transient heat conduction in polar coordinates with multiple layers in radial direction[J]. International Journal of Thermal Sciences，2008，47(3)： 261-273.

[6] Jain P K，Singh S，Rizwan-uddin. Analytical solution to transient asymmetric heat conduction in a multilayer annulus[J]. Journal of Heat Transfer，2009，131(1)： 73-92.

[7] 張國柱，夏才初，殷卓. 寒區隧道軸向及徑向溫度分布理論解[J]. 同濟大學學報(自然科學版), 2010, 38(8): 1117-1122. (ZHANG Guozhu, XIA Caichu, YIN Zhuo. Analytical solution to axial and radial temperature of tunnel in cold region[J]. Journal of Tongji University (Natural Science), 2010, 38(8): 1117-1122. (in Chinese))

[8] Battaglia J L, Toutain J, Pradere C, et al. Numerical inversion of Laplace transform for time resolved thermal characterization experiment[J]. Journal of Heat Transfer, 2011, 133 (4): 626-638.

[9] 邵珠山，喬汝佳，王新宇. 高地溫隧道溫度與熱應力場的彈性理論解[J]. 巖土力學, 2013(增刊1): 1-8. (SHAO Zhushan，QIAO Rujia，WANG Xinyu. Elasticity solution for temperature and stress fields of tunnels with high geothermal temperature[J]. Rock and Soil Mechanics, 2013(S1): 1-8. (in Chinese))

[10] 鄭陽，王道遠，朱永全. 寒區運營隧道溫度場數值模擬[J]. 公路工程, 2014, 39(1): 40-43. (ZHENG Yang, WANG Daoyuan, ZHU Yongquan. Numerical simulation of temperature field in cold region operational tunnel[J]. Highway Engineering, 2014, 39(1): 40-43. (in Chinese))

[11] Okada K. Prevention by adiatic treatment of tunnel lining[J]. Japanese Railway Engineering, 1985, 26 (3): 75-80.

[12] 陳建勛. 寒冷地區隧道防凍隔溫層設計計算方法及應用[J]. 土木工程學報, 2004, 37(11): 85-88. (CHEN Jianxun. Design method and application for anti-freezing thermal-protective layers in cold-area tunnels [J]. China Civil Engineering Journal, 2004, 37(11): 85-88. (in Chinese))

[13] 夏才初，范東方，韓常領. 寒區隧道不同類型凍土段隔熱(保溫)層鋪設厚度計算方法[J]. 中國公路學報, 2013, 26(5): 131-139. (XIA Caichu, FAN Dongfang, HAN Changling. Piecewise calculation method for insulation layer thickness in cold region tunnels[J]. China Journal of Highway and Transport, 2013, 26(5): 131-139. (in Chinese))

[14] 陳建勛，羅彥斌. 寒冷地區隧道溫度場的變化規律[J]. 交通運輸工程學報, 2008, 8(2): 44-48. (CHEN Jianxun, LUO Yanbin. Changing rules of temperature field for tunnel in cold area [J]. Journal of Traffic and Transportation Engineering，2008，8(2)：44-48. (in Chinese))

[15] 郝飛，孫全勝. 寒區公路隧道溫度場及保溫層的研究[J]. 現代隧道技術, 2012, 49(1): 39-43. (HAO Fei, SUN Quansheng. Research on the temperature fields and insulating layers of highway tunnels in a cold region [J]. Modern Tunnelling Technology, 2012，49(1): 39-43. (in Chinese))

[16] 周小涵. 通風對寒區隧道抗防凍的影響及合理設防范圍研究[D]. 成都: 西南交通大學, 2012: 29-31. (ZHOU Xiaohan. Study of the influence of ventilation on frost resistance and reasonable range of resistance for tunnel in cold region[D]. Chengdu: Southwest Jiaotong University, 2012: 29-31. (in Chinese))

[17] ZHOU Xiaohan, ZENG Yanhua, FAN Lei. Temperature field analysis of a cold-region railway tunnel considering mechanical and train-induced ventilation effects [J]. Applied Thermal Engineering, 2016, 100: 114-124.

[18] 張建榮，劉照球. 混凝土對流換熱系數的風洞實驗研究[J]. 土木工程學報, 2006, 39(9): 39-42. (ZHANG Jianrong, LIU Zhaoqiu. A study of the convective heat transfer coefficient of concrete in wind tunnel experiment [J]. China Civil Engineering Journal, 2006，39(9): 39-42. (in Chinese))

Finite Difference Calculation Model for Tunnel Temperature Field Based on Circular Cross-section

ZHOU Xiaohan1， 2, ZENG Yanhua1,*, BAI Yun1, WEI Yingjie3, HE Xing4

(1．KeyLaboratoryofTransportationTunnelEngineeringofMinistryofEducation,SouthwestJiaotongUniversity,Chengdu610031,Sichuan,China; 2．ChinaRailwayMajorBridgeReconnaissance&DesignInstituteCo.,Ltd.,Wuhan430050,Hubei,China; 3．ChengduConstructionEngineeringGroupCorporation，Chengdu610014,Sichuan,China; 4．ChinaRaiwayEryuanEngineeringGroupCo.,Ltd.,Chengdu610031,Sichuan,China)

The horseshoe-shaped tunnels and end wall type tunnels are common in mountain areas. It’s necessary to study the suitability of tunnel temperature field calculation model using circular cross-section. The temperature fields of cold region horseshoe-shaped tunnels in Northeast China are calculated by circular cross-section model using finite difference method of convection-conduction equations based on air-lining-surrounding rock. The calculation results are compared with measured data. The results show that: 1) The finite difference calculation method by using circular cross-section model and considering effects of wind flow velocity and wind flow temperature, is easy and convenient; it can meet the construction requirements. 2) The wind flow velocity in tunnel and wind flow field at tunnel portal have a significant influence on tunnel temperature field. The temperature of secondary lining surface increases by 7.2 ℃ when wind flow temperature at tunnel portal increases by 10 ℃ . The temperature of secondary lining surface reduces by 6.6 ℃, 2.7 ℃, 1.5 ℃ and 0.9 ℃ when the wind flow velocity is 2 m/s, 3 m/s, 4 m/s and 5 m/s respectively.

horseshoe-shaped tunnel; circular tunnel; temperature field; finite difference; field monitoring; wind flow temperature; wind flow velocity

2016-02-03;

2016-10-31

國家自然科學基金資助項目(51278426)

周小涵(1988—)，男，重慶忠縣人，西南交通大學橋梁與隧道工程專業在讀博士，從事大型交通隧道的防災救援技術和地鐵結構研究。E-mail: zhouxh2008@126.com。 *通訊作者： 曾艷華， E-mail: zengyhua@163.com。

10.3973/j.issn.1672-741X.2016.11.007

U 45

A

1672-741X(2016)11-1332-05