, 张泽, 高钰 东北大学 资源与土木工程学院, 辽宁 沈阳 110819
收稿日期：2018-07-30
基金项目：国家自然科学基金资助项目(51474050，U1602232)；地质灾害防治与地质环境保护国家重点实验室项目(SKLGP2014K011)；辽宁省高等学校优秀人才支持计划项目(LN2014006)。
作者简介：杨天娇(1993-)，女，辽宁锦州人，东北大学博士研究生;
王述红(1969-)，男，江苏泰州人，东北大学教授，博士生导师。
摘要：利用“三区域”理论, 建立考虑水分迁移和水冰相变的寒区隧道水热耦合问题的联合求解微分方程, 然后采用COMSOL Multi-physics软件实现温度场和水分场耦合数值模拟, 进而将数值模拟结果与土柱冻结试验的结果进行对比, 验证水热耦合数值模拟模型的有效性.根据冻胀理论, 着重对寒区隧道的应力场控制方程进行推导, 利用孔隙冰与冻胀率之间的关系建立寒区隧道水热力三场耦合计算模型.最后以牡丹江到绥芬河区间的牡绥隧道为例, 对温度场、水分场以及应力场进行了模拟分析.结果表明:外界温度的变化对开挖后的隧道会有很大影响.随着时间的增加, 隧道洞口冻结圈厚度逐渐增加, 在1月份达到最大冻深2 m左右.
关键词：寒区隧道水热耦合水热力耦合冻胀COMSOL
Numerical Analysis on Thermo-Hydro-Mechanical Coupling of Surrounding Rocks in Cold Region Tunnels
YANG Tian-jiao, WANG Shu-hong
, ZHANG Ze, GAO Yu School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China
Corresponding author: WANG Shu-hong, E-mail: shwangneu@126.com
Abstract: Based on the theory of "three regions", a combined differential equation is established to solve the problems of hydro-thermal coupling in cold regions considering water transfer and water ice phase transition. The numerical simulation of temperature field coupled with water field was carried out by COMSOL multi-physics software, and the numerical simulation results were compared with the results of soil column freezing experiment to verify the effectiveness of the hydrothermal coupled numerical simulation model. By considering the theory of frost heave, the stress field governing equation of the tunnel in the cold regions is deduced, and the thermo-hydro-mechanical coupled model of the tunnel in the cold regions is established by using the relationship between pore ice and frost heave rate. Finally, the temperature field, humidity field and stress field were simulated by taking the Musui tunnel between Mudanjiang and Suifenhe as a case study. The results show that the change in external temperature has a great impact on the tunnel after excavation. The thickness of the freezing ring of the tunnel hole gradually increased with time and reached a maximum depth of about 2 meters in January.
Key words: cold region tunnelhydro-thermal couplingthermo-hydro-mechanical couplingfrost heaveCOMSOL
我国有一半以上的国土属于寒区, 而且大部分分布在东北和西部地区[1], 所以这些地方的工程建设常常会受到低温环境的影响, 冻害现象非常严重.寒区隧道修建过程中当温度降到0℃以下时围岩会受到冻胀的影响, 存在于岩体中的原位孔隙水及来自远处的补给水都会由于温度低于冰点冻结成冰, 从而造成岩土体的体积膨胀, 在隧道的结构上产生冻胀力.这个过程是一个复杂的水热力耦合问题.长期以来, 不少学者采用水热耦合模型针对冻土区和未冻区两个区域研究寒区隧道冻胀问题[2-3].Konrad和Nakano等[4-5]分别从实验和理论方面证实在未冻土和冻土区之间存在一个冻结缘或正冻区.谭贤君等[6]对Nakano提出的“三区域”理论进行了饱和岩土体的数值分析, 指出他并没有考虑岩土体中水分迁移对温度场的影响和冰对未冻水的阻碍作用.Neaupane和Gatmiri等[7-8]建立了低温冻结岩体THM耦合控制方程, 并通过编制有限元程序进行模拟分析.
本文首先基于三区域理论, 对原有模型进行改进并应用于非饱和土的低温水热耦合问题, 考虑水分迁移对寒区温度场的影响，并引入阻抗系数以考虑冰对未冻水的阻碍作用.使用COMSOL软件[9]进行二次开发并与文献[10]中经典试验进行对比验证, 证明本文所建立的水热耦合数学模型的正确性.然后基于冻胀理论, 着重对寒区隧道的应力场控制方程进行推导，建立寒区隧道水热力三场耦合计算模型, 从而对温度场、水分场以及位移场进行模拟分析, 最后对隧道围岩的导热系数、体积热容和渗透系数进行参数敏感性分析, 研究成果能较真实地反映寒区隧道冻害现象发生过程, 具有一定的参考价值.
1 寒区隧道水热耦合以现有的热平衡原理Harlar模型为基础, 结合岩土介质实际冻结过程, 基于三区域理论, 考虑水冰相变和水分迁移的影响, 提出了三区域分析模型, 建立了更符合实际、更全面的寒区隧道水热耦合模型.
1.1 基本假设1) 土介质为均质各向同性孔隙介质, 由未冻水、冰、岩体骨架、空气组成, 冻结过程中不发生变形;
2) 岩土介质的冻结和融化过程中忽略水气转化和迁移, 忽略空气对水分迁移的贡献;
3) 忽略水分流动过程的对流传热和水分蒸发, 只考虑水冰相变的潜热和热传导过程.
1.2 温度场控制方程根据傅里叶定律和能量守恒定律, 分别在未冻区、已冻区和相变域建立考虑热传导和对流扩散项的微分控制方程.
已冻区:
 | (1) |
 | (2) |
 | (3) |
 | (4) |
 | (5) |
为了方便计算, 把三区域控制方程统一写成一个方程式的形式, 通过调节参数来达到在三区域之间的转化.
 | (6) |
低温下土中的水有两个冻结过程:原位水冻结和水分迁移冻结.令冰水相变点的温度为θm, 则
 | (7) |
 | (8) |
 | (9) |
在冻土中, 孔隙冰的存在占据了岩土体中部分的孔隙体积, 使得能够导水的孔隙减少, 导水通道的减少削弱了未冻水分的迁移.为了体现与常温下非饱和土水分迁移的区别, 考虑了负温下孔隙冰的阻抗作用, 引入阻抗系数I，I=1010θi.联立式(6)和式(9)得到考虑水分迁移的寒区隧道围岩的水热耦合方程.
2 模型验证为了验证数学模型的正确性和可行性, 应用文献[10]中的试验进行反演验证.试验为封闭系统下一维垂直土柱冻结试验, 如图 1所示.土样为张掖壤土, 由非饱和土组成, 采用从上至下的无外载垂直冻结方式, 试验参数为:高度13.68 cm, 干密度1 500 kg/m3, 初始含水量0.220 8, 孔隙比0.444.
图 1(Fig. 1)
 | 图 1 土柱试验装置示意图Fig.1 Diagram of soil column test device |
在文献[10]试验中, 并没有给出热学参数, 所以本文热学参数参照文献[11]中粉土的热参数.未冻区和已冻区土壤的导热系数分别为λu=1.13W/(m·℃), λf=1.58W/(m·℃)；未冻区和已冻区的比热容分别为cu=2.36 MJ·m-3·℃, cf=2.82 MJ·m-3·℃.
根据文献[10]中的实测数据分别得到在冻土和未冻土中的扩散系数D和渗透系数K随饱和度S变化的函数如下:
 | (10a) |
 | (10b) |
表 1(Table 1)
 表 1 试验初始温度Table 1 Experimental initial temperature
| 表 1 试验初始温度 Table 1 Experimental initial temperature |
图 2(Fig. 2)
 | 图 2 土柱上下边界温度变化Fig.2 Temperature variation of the upper and lower boundary of the soil column |
分别计算0.6, 8.8和47.2 h的温度场变化情况, 并与文献[10]中的试验和数值模拟进行了对比, 结果如图 3和图 4所示, 基本吻合较好, 验证了模型的正确性.
图 3(Fig. 3)
 | 图 3 温度场变化对比图Fig.3 Contrast diagram of temperature field change |
图 4(Fig. 4)
 | 图 4 冻结深度变化对比图Fig.4 Contrast diagram of frozen depth change |
3 寒区隧道水热力耦合3.1 水热力耦合方程冻胀的产生是由于岩土介质中孔隙水在低温状态下发生相变, 水变成冰导致的体积增大.这种体积膨胀与材料的热膨胀现象相似, 所以本文利用热膨胀的计算理论来解释此冻胀模型.
热膨胀引起的应变为εi=α(θ－θ0).式中:α为材料的线膨胀系数; θ为瞬时温度(℃), θ0为初始温度(℃).
寒区隧道围岩的冻胀一般用冻胀率来反映, 它等于围岩冻结后膨胀的体积ΔV与围岩冻结前的体积V之比, 实际上在侧限条件下就是前后高度的变化.
体积冻胀率为
 | (11) |
 | (12) |
 | (13) |
3.2 水热力耦合模型将以上水热力耦合数值仿真方法应用于牡丹江到绥芬河区间的牡绥隧道分析.选取隧道未开挖前的上覆山体和隧道下一部分山体进行计算.隧道牡丹江端洞口属缓坡地貌, 地形坡角在5°~15°；隧道绥芬河端属缓坡地貌, 地形坡角在5°~10°.计算模型如图 5所示.
图 5(Fig. 5)
 | 图 5 山体计算模型Fig.5 Calculation model of mountain |
3.2.1 牡绥隧道初始温度场山体的初始温度按0.04 ℃/m的梯度计算.山体表面受环境温度影响较大, 所以山体上表面为温度边界, 具体形式根据当地年温度变化进行拟合.
根据现场测试, 得到山体的导热系数和比热容随温度变化的曲线，如图 6所示.
图 6(Fig. 6)
 | 图 6 导热系数和比热容的变化Fig.6 Variation of thermal conductivity and specific heat capacity |
黑龙江省牡丹江市一年的温度变化情况如图 7所示, 大致认为一年内的大气温度是随着时间呈现周期性变化的, 一般写成正弦或余弦函数的形式, 如式(14)所示.
 | (14) |
 | 图 7 牡丹江市一年温度变化图Fig.7 Temperature chart of one year in Mudanjiang |
式中:θ为温度(℃); t为时间(s).
分别对隧道区最冷月份和最热月份隧道未开挖前的山体温度场进行数值模拟, 模拟结果如图 8所示.从图中可以发现不管是最冷的1月还是最热的7月, 隧道轴线附近的温度场变化都不大, 绝大部分为恒温场, 温度场不随季节的变化而变化, 只有在山体表面一定厚度内, 温度场发生改变.
图 8(Fig. 8)
 | 图 8 温度场Fig.8 Temperature fields (a)—7月；(b)—1月. |
从计算结果可以看出，隧道未开挖之前在隧道轴线位置的温度场基本不受外界环境的影响, 但是随着隧道的开挖, 隧道洞口受外界寒冷气温的影响, 往往会发生冻害现象, 所以着重对隧道洞口的温度场、水分场和应力场进行计算分析.
3.2.2 牡绥隧道洞口耦合计算模型取隧道牡丹江端洞口为主要研究对象, 进行耦合计算与分析, 洞口计算模型如图 9所示.
图 9(Fig. 9)
 | 图 9 隧道洞口模型示意图Fig.9 Diagram of tunnel portal model |
牡丹江端洞口围岩主要为粉质黏土并在洞口喷射60 cm厚的C25混凝土衬砌, 围岩及衬砌的参数见表 2, 其中n为孔隙率，θf为洞口温度.
表 2(Table 2)
表 2 模型计算物理参数Table 2 Model calculation physical parameters
| 表 2 模型计算物理参数 Table 2 Model calculation physical parameters |
围岩的导热系数和体积热容是随温度变化的函数, 具体变化形式如表 3所示.
表 3(Table 3)
表 3 模型热学参数Table 3 Model thermal parameters
| 表 3 模型热学参数 Table 3 Model thermal parameters |
根据对隧道所在山体的初始温度场计算, 可知隧道洞口初始温度为12 ℃.洞口的温度根据外界环境温度的变化而变化, 洞口温度边界参考牡丹江市全年温度变化曲线, 见图 7.洞口围岩初始饱和度为0.8.
3.3 牡绥隧道洞口水热力耦合模型分析3.3.1 温度场分析根据隧道所在地的年变化温度, 计算一年内隧道洞口处围岩的温度场变化情况.本文重点给出了11月、12月、1月和2月的围岩温度场以及隧道洞口处冻结圈的变化情况.
从图 10中可以看出，从11月开始, 随着时间的增加, 冻结圈的厚度也在增加, 在1月达到最大, 最大冻深为2 m.从图 11a中也可以看出最冷温度出现在1月, 在2月的温度场分布图中可以明显看出有温度回暖现象.
图 10(Fig. 10)
 | 图 10 11月到来年2月的冻结圈Fig.10 Frozen circles from November to February next year |
图 11(Fig. 11)
 | 图 11 11月到次年2月水热力耦合模型分析Fig.11 Analysis of thermo-hydro-mechanical coupling model from November to February next year (a)—温度场分布；(b)—含冰量分布；(c)—冻胀力分布. |
3.3.2 水分场分析在东北地区, 11月、12月、1月、2月是低温月份, 室外温度处于零度以下, 为水冰相变提供了条件, 也是隧道冻胀现象的高发月份, 所以着重对这4个月进行分析.
从图 11b中可以看出, 随着季节和大气温度的变化, 隧道洞口围岩的含冰量也在变化, 1月时含冰量达到峰值.
3.3.3 应力场分析求解出11月—次年2月隧道洞口含冰量后, 利用所建立的水热力三场耦合模型对隧道洞口这四个月的冻胀力进行计算, 结果如图 11c.可以看出最冷的1月的冻胀力达到最大值，为1.4 MPa, 2月时温度逐渐回升, 冻胀力也有所减小，为1 MPa.
4 结论1) 在原有水热耦合模型基础上, 推导出考虑水冰相变的三区域温度场方程和水分迁移的水分场方程, 从而建立了更符合实际的水热耦合模型.利用COMSOL软件与一维垂直土柱冻试验进行验证, 数值模拟得到的温度场和冻结深度的图像与试验结果吻合较好, 进而验证了模型的正确性.
2) 在建立的水热耦合模型的基础上, 考虑冻胀效应的影响, 利用孔隙冰与冻胀率之间的关系建立了寒区隧道水热力耦合的计算模型.以COMSOL Multiphysics多场耦合软件为平台, 求出隧道所在山体的初始温度场, 结果表明, 主要是开挖通风后洞口受到环境温度影响较大, 常常会发生冻害现象.
3) 本文建立的水热力三场耦合模型对牡绥隧道牡丹江端隧道洞口11月、12月、1月和2月的温度场、水分场和应力场进行计算分析, 通过计算结果可以看出，随着时间的增加, 隧道洞口冻结圈厚度逐渐增加, 1月达到最大冻深2 m左右.含冰量峰值出现在1月, 同时冻胀力达到最大值1.4 MPa.
参考文献
| [1] | 景远吉.寒区隧道围岩水热力耦合数值分析[D].西安: 西安科技大学, 2017. ( Jing Yuan-ji.The moisture-heat-stress coupling numerical analysis of the tunnel surrounding rock in the cold regions.[D]. Xi'an: Xi'an University of Science and Technology, 2017.http://cdmd.cnki.com.cn/Article/CDMD-10704-1017725511.htm) |
| [2] | Harlan R L. Analysis of coupled heat-fluid transport in partially frozen soil[J].Water Resources Research, 1973, 9(5): 1314–1323.DOI:10.1029/WR009i005p01314 |
| [3] | Taylor G S, Luthin J N. A model for coupled heat and moisture transfer during soil freezing[J].Canadian Geotechnical Journal, 1978, 15(4): 548–555.DOI:10.1139/t78-058 |
| [4] | Konrad J M, Morgenstern N R. Effects of applied pressure on freezing soils[J].Canadian Geotechnical Journal, 1982, 19(4): 494–505.DOI:10.1139/t82-053 |
| [5] | Nakano Y. Quasi-steady problems in freezing soils(Ⅰ):analysis of the steady growth of an ice layer[J].Cold Regions Science and Technology, 1990, 17(3): 207–226.DOI:10.1016/S0165-232X(05)80002-3 |
| [6] | 谭贤君, 余祥宏, 陈卫忠, 等. 岩土介质在冻融过程中的温度场研究及工程应用[J].岩土工程学报, 2012, 31(sup1): 2867–2874. ( Tan Xian-jun, Yu Xiang-hong, Chen Wei-zhong, et al. Temperature field research and engineering application of geotechnical media during freezing and thawing process[J].Chinese Journal of Geotechnical Engineering, 2012, 31(sup1): 2867–2874.) |
| [7] | Neaupane K M, Yamabe T. A fully coupled thermo-hydro-mechanical nonlinear model for a frozen medium[J].Computers and Geotechnics, 2001, 28(8): 613–637.DOI:10.1016/S0266-352X(01)00015-5 |
| [8] | Gatmiri B, Delage P A. Formulation of fully coupled thermal-hydraulic-mechanical behaviour of saturated porous media-numerical approach[J].International Journal for Numerical & Analytical Methods in Geomechanics, 2015, 21(3): 199–225. |
| [9] | COMSOL.COMSOL multiphysics user's guide[R]. Version 4.3b.Stockholm: COMSOL AB, 2012. |
| [10] | 胡和平, 杨诗秀, 雷志栋. 土壤冻结时水热迁移规律的数值模拟[J].水利学报, 1992(7): 1–8. ( Hu He-ping, Yang Shi-xiu, Lei Zhi-dong. Numerical simulation of hydrothermal migration in soil freezing[J].Journal of Hydraulic Engineering, 1992(7): 1–8.DOI:10.3321/j.issn:0559-9350.1992.07.001) |
| [11] | 徐学祖, 王家澄, 张立新. 冻土物理学[M]. 北京: 科学出版社, 2010. ( Xu Xue-zu, Wang Jia-cheng, Zhang Li-xin. Frozen soil physics[M]. Beijing: Science Press, 2010.) |
