寒区隧道围岩水热耦合数值分析

张 泽， 王述红， 杨天娇， 张雨浓

(东北大学 资源与土木工程学院， 辽宁 沈阳 110819)

我国东北部以及西部这些高纬度或高海拔地区有一半以上的国土属于寒区[1].寒区隧道在修建过程中由于温度降低围岩会受到冻胀的影响，除了岩土体孔(裂)隙中的原位水发生冻结外，岩土体中的水分迁移使得水分不断向冻结锋面迁移，促使冻胀更加明显.因此，隧道岩体中的冻胀问题归根结底为土体中的温度水分迁移问题.冻融状态转变过程中水分迁移导致土体中水分占比变化引起土体热参数发生改变；融化界面的冰晶会对周围的水分产生吸引作用从而导致界面水头的产生；由于水冰相变土体中水分含量降低而冰含量相应升高因此阻碍水的扩散，这些会对冻结过程造成深远影响；同时相变释放的热量也会对冻结产生较大影响；即使在冻结温度时土体中仍然存在一部分未冻水，且其含量与负温保持动态平衡关系.

隧道冻害过程是温度-渗流耦合问题，长期以来，各国学者建立了不同的温度-渗流耦合模型对寒区隧道的冻融进行研究[2-4].Liu等[5]采用热力学、流体力学、弹塑性力学基本理论，建立了非饱和多孔介质材料低温THM耦合模型，对冻胀量的分布、冻胀深度以及应力场、温度场的分布进行模拟计算分析.谭贤君等[6]通过数值模拟对提出的“三区域”理论进行分析，但是对于水分迁移以及水冰相变的阻碍作用没有做出解释.Yamabe等[7-8]对饱和的砂岩和凝灰岩进行了实验并对低温水冰相变下的饱和岩石多次耦合过程进行了研究.徐光苗等[9]首次对低温岩石进行研究，建立了低温岩石的THM三场耦合方程，并采用有限元程序进行模拟.

本文基于混合物理论，考虑水分迁移以及水冰相变的影响，以及土体孔隙冰对未冻水的阻滞作用对原有模型进行改进并应用于非饱和土的低温水热耦合问题中，并使用COMSOL软件[10]进行模块开发与文献[11]中经典实验进行数值验证，证明本文所建立的温度-渗流耦合数学模型的适用性.然后以西藏自治区米林隧道为工程实例，对考虑水分迁移与不考虑水分迁移的温度场、水分场进行了模拟分析，研究成果能较真实地反映寒区富水隧道冻害现象发生过程，具有一定的参考价值.

1 寒区隧道水热耦合

为了建立更符合寒区隧道冻胀问题的温度-渗流耦合模型，在现有经典Harlar水热模型基础上，将水分迁移和水冰相变考虑其中建立更为全面的岩土介质冻结过程的温度-渗流耦合模型，为寒区隧道的冻害问题提供参考.

1.1 基本假设

1) 土介质为均质各向同性孔隙介质，由未冻水、冰、岩体骨架组成，冻结过程中不发生变形.

2) 岩土介质中水分迁移，忽略空气对水分迁移的贡献.

3) 只考虑水冰相变的潜热和热传导过程.

4) 不考虑由于冰透镜体等冰体形成造成的上部已冻土克服下部水分重力对水分的抽吸作用.

1.2 温度场控制方程

根据傅里叶定律和能量守恒定律，在温度场控制方程中的比热容、导热系数都为含水率的函数，结合冻土的相变条件，单元体内的体积热容可描述为单元体的冰水相变所产生或吸收的相变热，由此可得

(1)

从而可得温度场的控制方程：

(2)

式中：Ca为等效体积热容量，kJ/(m3·℃)；t为土体瞬态温度,℃；τ为时间,s；λ为热传导系数[W/(m·℃)]；ρw,ρi分别为水、冰密度,kg/m3；θi为冰体积含量；L为相变释放潜热，取值334 kJ/kg；为哈密顿算子.

1.3 水分场控制方程

在寒区隧道工程围岩中即使在冻结点温度以下仍然存在没有冻结的水分，其遵循达西定律进行迁移：在单元体内考虑孔隙冰对未冻水的阻滞作用，非饱和土中的未冻水迁移微分方程即为水分场耦合温度控制方程：

(3)

式中,D(θu)，K(θu)分别为非饱和土扩散系数和导水率，可见其都是随着未冻水含量发生变化.

冻土中水的扩散率D(θ)计算：

(4)

式中：k(θ)为非饱和土的渗透率,m/s；c(θ)为比水容量,1/m，可由滞水模型确定.

在冻土中，孔隙冰含量和孔隙水含量存在动态变化，冰的增加必占据孔隙从而减少水分迁移的通道.为了体现与常温下非饱和土水分迁移的区别，考虑了负温下孔隙冰的阻抗作用，引入阻抗系数的概念.I为阻抗系数.

I=10-10θi，

(5)

k(θ)=ksSl(1-(1-S1/m)m)2,

(6)

c(θ)=a0m/(1-m)·S1/m·(1-S1/m)m.

(7)

其中：l为本构参数；S为冻土相对饱和度，定义为

(8)

式中,θr和θs分别代表土的剩余含水量和饱和含水量.

2 关键耦合参数描述

1) 等效热容:

(9)

2) 等效热传导系数.导热系数的取值有多种加权方法，本文采用指数加权方法，可表示为

λ(θ)=(λs)θs(λl)θl(λi)θi.

(10)

3) 未冻水含量.为真正建立水分场与温度场的耦合方程，需要引入θ,θi与t之间的联系方程.徐学祖等[12]基于试验数据，给出了冻土中未冻水含量与温度的经验关系式：

(11)

式中：tf为土体冻结温度,℃；w0为土体初始含水量,%；wu为负温度为t时的未冻水质量分数,%；B是与土类和含盐量确定的常量.按照徐学祖等[12]冻土中水分迁移的实验研究中的测定方法，砂土、粉土、黏土中的B值为0.61，0.47，0.56.

4) 固液比.冻土中孔隙冰体积与未冻水体积之比，记为Bi:

(12)

因此，冻土中孔隙冰、未冻水和温度的联系方程为

θi=Bi(t)·θu.

(13)

3 水热力模型的试验验证

文献[11]中的一维土柱冻结试验可以对所建立的模型正确性和可行性进行验证，该试验为封闭系统,如图1所示.

土样由非饱和土组成，无外载冻结方式为从上至下，试验参数见表1.

表1 试验参数

文献[11]中给出了试验的初始温度见图2，土柱上边界和下边界温度的变化值见表2.

表2 试验上下边界温度变化

本文计算的物理参数如表3所示.

表3 计算物理参数

根据文献中的实验结果分别对528 min和2 830 min 的温度变化进行计算，与文献[11]中的试验数值进行了对比，如图3所示，验证结果基本吻合，证明了模型的正确性.

4 工程实例

4.1 工程概况及计算参数

米林隧道址山势雄伟，区内高点位于隧道轴线左侧山脉，海拔高度为4 230 m，最低点位于隧道出口宽谷地带，标高为2 940 m，典型的高海拔隧道，最大埋深为1 200 m.项目所在地区历年年平均气温4.3 ℃，历年极端最高气温为21 ℃，历年极端最低气温为-25 ℃，历年最冷月平均气温为-7.6 ℃，按对隧道工程影响的气候分区，属寒冷地区.米林隧道为高海拔富水寒区大长隧道.

根据地质勘查资料，选取隧道的上覆山体和隧道下一部分山体进行计算，取其硐深10 m处的断面为研究对象，该断面处的地质条件为：地质成分主要为块石土、粉土.根据隧道的埋深以及断面尺寸，建立隧道均匀孔隙介质简化计算模型，如图4所示.

1) 米林隧道温度场初始条件和边界条件.

林芝县近五年月平均气温的变化情况：林芝县全年的平均气温为2.3 ℃，最冷月为一月,平均气温为-9.8 ℃，最热月为7月,平均气温为11.8 ℃；极端最高气温为27 ℃，极端最低气温为-13 ℃，由于海拔高，温度垂直变化明显，海拔每升高100 m温度下降0.6 ℃.

将五年的每月平均温度绘制如图5所示,大致认为一年内的温度变化规律大致呈现周期性变化，可写成正弦或余弦函数的形式:

(14)

实际中隧道顶部和隧道内部的大气温度变化会引起围岩中的温度-渗流变化.因此选取隧道断面施工期的地面(ab边)及内部(efg边)温度为大气温度变化，见式(14).根据对称性ae,gd边为绝热边界，同样bc,cd边也为绝热边界.由于为高海拔隧道，需要考虑高度对于温度的影响，隧道顶部温度要低于隧道围岩温度如图6所示.根据施工现场提供的隧道施工阶段温度值，隧道围岩初始温度为-0.98 ℃.

本文主要研究低温状态下隧道中水热耦合效应，因此选取本年11月到次年4月这六个月进行研究，温度取值变化趋势如图6所示.

2) 渗流场初始边界条件与边界条件.

bc边为水头压力边界，若将d点作为原点进行计算，水压力在竖直方向上呈线性变化，总水头等于基质吸力水头与位置水头之和，因此取压力水头为式(15)所示，初始含水量根据现场提供数据取值为12%.

(15)

式中：p为水头压力，Pa；h为高度，m.

根据对称性ae，gd为不透水边界，衬砌与围岩的接触边efg压力水头取值为0，cd边为固定边界，压力为0，忽略衬砌的渗透性，初始水压力取值为0.

由地质条件可知，米林隧道计算所取断面主要为粉土，并喷射60 cm厚的C25混凝土衬砌.主要计算物理参数如表4所示，其他参数如表3所示.

表4 模型计算主要参数

4.2 米林隧道水热耦合模型分析

1) 温度场时空分布.林芝地区在11月份开始进入零下温度，并持续降温.图7显示了隧道整体由11月份到次年4月份经历的温度变化状况.根据当地的气温变化规律，隧道边界气温持续降低，隧道顶部边界气温由-0.82 ℃降低到-9 ℃，隧道内部边界温度由-0.74 ℃降低到-11.11 ℃，导致隧道整体温度的下降，下降趋势为由顶部向底部，由洞口向内部延伸扩展.

距离边界越近的岩体温度下降越快，这是因为边界带范围内的温度梯度大，在导热系数相同的情况下，热量传递量大；此外再考虑冻结状态下水冰相变的情况， 水相变为冰时比热容会降低而导热系数会增大，由此也增加热量的传输速率.在3月和4月温度有明显的回升.

2) 未冻水含量时空分布.由上述温度时空分布的分析知，在11月到次年2月的120 d中隧道顶部及隧道内部经历降温的过程，隧道岩体存在局部降温，降温范围可以从图8a～图8d观察得出，红色区域持续减小，说明此区域存在持续降温状态.

由图8e,图8f分析可得，由于隧道顶部及内部温度升高导致低温时已经冻结成冰的固态水转化为液态水， 显示了隧道顶部边界和隧道边界由模拟三月时逐渐转变为正温，说明边界进入了融化状态，导致隧道顶部边界和隧道内部边界未冻水含量升高.在隧道顶部和隧道内部未冻水含量由120 d时的5.3%迅速增加到180 d时的9.8%，这对隧道的稳定性是十分不利的.

3) 冰含量时空分布.由图9a～图9d可以看出，隧道岩体的固态冰含量持续升高，这是因为隧道岩体持续降温，液态水被冻结为固态水，从而升高了固态冰的含量，由上述温度时空分布的分析知，在11月至2月隧道顶部以及隧道内部经历降温过程，隧道岩体中存在局部降温，降温的范围可以从图9a～图9d观察出来，蓝色区域(固态冰含量5%～10%)持续减小，由蓝绿色区域的扩张代替(固态冰含量10%～25%)，说明此区域经历着持续降温状态，并且含冰量峰值出现在1月份.

由图9e，图9f分析可得，由于隧道顶部以及隧道内部的升温导致已经冻结成冰的固态水转为液态水，因此隧道在3月份开始内部的冻结冰开始融化，在4月份时随着温度的持续升高，其含冰量也随之降低，由冻结时的25%降低到15%，可以推测随着温度的持续升高，隧道岩体中的冰含量将一直减小直到全部融化为液态水.岩体中的水冰相变和冰体的冻结融化对于隧道整体的稳定性产生重要的影响.

4) 水分迁移对温度场的影响分析.为与上述已经考虑水分迁移的模型对比，应用式(16)建立不考虑水分迁移的数值模型，只考虑隧道岩体整体的热迁移，不考虑冰水相变.

(16)

通过对比图10和图11可见，未考虑水分迁移的温度场中热传导速度较快，图10的红色区域较图11的红色区域减小的速率大，其黄绿色区域(-4 ℃～-7 ℃)扩张速度大，这是因为考虑水分迁移的热传导过程受到的水相变成冰释放潜热的影响，相变潜热提供的热量阻止了低温的传递，而且水的比热容一般是岩体的比热容的4倍左右，说明岩石升高同等温度要消耗4倍的热量，也在一定程度上阻碍了隧道岩体中的温度升高过程.

由3月份温度场对比分析可以看出，考虑水分迁移的隧道顶部以及隧道内部温度场与不考虑水分迁移的温度场相比，前者升温过程也比较缓慢，分析原因，导致这种情况主要有两点：一是由于水的比热容较大，升高相同温度消耗的热量较岩体的多；二是随着温度周期性变化，由于岩体升温导致岩体内的冰含量减少，而冰融化成水是吸收热量的过程，要消耗热能，相对的水成冰要放热，冰成水要吸热，所以考虑水分迁移的隧道岩体在升温过程中的传热较慢.

5 结 论

1) 在原有水热耦合模型基础上，推导出考虑水冰相变和水分迁移的水热耦合方程，从而建立了更全面更符合实际的水热耦合模型.利用COMSOL软件开发新的模块与一维土柱冻结实验进行反演验证，数值模拟得到的温度场与实验结果基本吻合较好，进而验证了模型的正确性.

2) 利用水热耦合模型对米林隧道11月到次年4月的温度场和水分场进行模拟分析.结果表明：隧道顶部边界气温由-0.82 ℃降低到-9 ℃，隧道内部边界温度由-0.74 ℃降低到-11.11 ℃，并在3月份隧道温度回升；同时隧道中冰含量在1月份达到峰值，在3月份开始减小，表明在隧道中存在明显的冻融现象，极易产生冻害问题.

3) 对比分析水分迁移对温度场的影响，考虑水分迁移的隧道顶部以及隧道内部温度场升温过程比较缓慢，证明固态冰的形成或融解时所产生的相变潜热对隧道中温度场的分布影响远远大于液态水靠重力迁移造成的热对流传热.