基于圆形断面的隧道温度场有限差分计算模型

周小涵， 曾艳华， 白 赟， 魏英杰, 何 省

(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