基于对流-导热耦合模型的寒区水工隧洞围岩温度-应力耦合分析

2021-11-12姜海波

水力发电 2021年8期

孟尧，姜海波

(石河子大学水利建筑工程学院，新疆石河子 832000)

0 引言

高海拔寒区隧洞在自然通风条件下寒冷空气与隧洞围岩发生对流换热，围岩之间也存在着导热，从而产生明显的温度差且出现温度应力，使围岩的应力分布及稳定性发生变化。寒区隧洞的安全稳定离不开隧洞温度-应力耦合的研究分析，因此有必要基于对流-导热耦合模型，对寒区水工隧洞围岩在不同环境条件下围岩不同位置的温度-应力耦合进行分析。目前，关于寒区隧洞围岩的温度场及温度-应力耦合的研究已取得很多成果。丁浩、孟尧等[1-2]研究发现，洞内外大气和围岩发生热传递，得到隧洞径向和纵向温度场的变化规律。严健等[3]研究得出，升温对掌子面周围的洞壁温度影响较大。徐衍、刘江川等[4-5]研究得出，围岩径向距离增加，温度影响逐渐消失，围岩的温度变化是影响围岩受力状态的重要因素。刘杰等[6]对比不同洞型开挖的应力场和塑性区分布，确定圆形断面为最优开挖断面。李松磊等[7]研究发现，隧洞在无衬砌时通水，温度荷载会产生较大的拉应力，影响工程的安全稳定。刘春龙等[8]研究得出，变温区温度场与隧洞围岩自重应力场的应力耦合关系曲线。王伟[9]研究发现，在对流荷载作用下，距离隧道围岩越近围岩温度变化越大，对流荷载可使隧道开挖后围岩应力减小，对隧道不同位置围岩应力的影响不同。王睢等[10]推导出深埋含衬砌有压圆形隧洞在施工和运行不同工况下的弹塑性解，验证了模型的正确性。贺耕夫等[11]研究发现，塑性区半径对围岩单轴抗压强度在一定范围内的变化非常敏感，甚至达到8～10倍洞径。汤盛显等[12]利用弹塑性交界面上应力连续原则，使用隧道壁面应变为边界条件，求出塑性区半径和对应的径向应变。王明年等[13]得出隧道洞周径向位移与强度比呈负相关，塑性区半径、隧道洞周径向位移均与开挖半径呈线性正相关的结论。刘思杰等[14]对开挖期间围岩破坏机理研究发现，在顶拱、拱肩多个部位出现应力集中区。任喜平[15]研究得出羊曲水电站泄洪洞围岩位移变形及应力分布规律。马康等[16]研究发现，随着温度荷载的增加，隧道产生径向应力与环向应力，径向应力在衬砌与围岩交界处出现应力集中现象。陈张华等[17]研究软岩力学参数的敏感性发现，降低软岩变形模量后厂区的位移、应力状态、塑性区均有变化。王平等[18]研究得出，深部高地应力软岩巷道塑性区轴比与应力集中是恶性循环的，提出控制巷道塑性区最大半径方向上的围岩变形是关键。李志雨[19]得出地下隧洞只在温度作用下以及在温度和原始地应力共同作用下隧洞围岩应力的数值解，并与解析解对比验证，两者结果吻合较好。以上研究分析了不同环境下寒区隧洞温度场的变化规律，并且进一步研究了温度荷载产生的温度应力作用下隧洞围岩发生的变化。

高海拔寒区隧洞贯通后在自然通风的低温条件下，围岩的力学性质发生变化，从而影响水工隧洞围岩的应力分布及稳定性。因此，有必要研究不同自然通风温度和不同风速下的水工隧洞洞口及洞中位置围岩温度-应力耦合及其演化规律。本文以新疆某寒区水工隧洞为依托，基于M-C本构模型，采用有限元仿真计算，对不同自然通风温度和不同风速下的水工隧洞洞口及洞中位置围岩温度-应力耦合及其演化规律进行系统分析，定量分析计算寒区隧洞温度-应力耦合的分布规律及其变化规律，研究隧洞内风温、风速对围岩温度-应力耦合的影响，并分析其影响机制，得到围岩瞬态的应力及塑性应变分布规律，可为寒区水工隧洞安全运行提供理论依据。

1 数值模型建立

1.1 工程实例

新疆某引水隧洞位于高寒半干旱气候区，季节性温差较大，最低气温-34.3 ℃，最高气温35.9 ℃，最冷月平均气温-16.66 ℃；多年平均风速2.1 m/s，风向西。根据现场监测得知，引水隧洞周围的围岩实测温度为9 ℃。洞址区为第四系地层，未见层理，洞室围岩完整性较好，上覆盖岩层厚260～300 m。通过地质勘查，隧洞区只有少量上层滞水，围岩及周围岩体水分较少，可忽略水对应力的影响。

1.2 对流-导热耦合模型构建及参数

计算流体与固体壁面接触换热用牛顿冷却公式表达，即

q=h(Ta-Tb)

(1)

式中，h为物体与周围介质的对流换热系数；Ta为物体的壁面温度；Tb为介质的平均气温。热传导用傅立叶导热基本定律表达

(2)

式中，q为单位传热面积热流量；n为等温面法线方向向量。对流-导热公式表示为

(3)

隧洞直径3 m，根据实际工程，采用有限元软件建立模型，围岩计算范围取为21 m×21 m×1 000 m的立方体模型。因为温度荷载对隧洞围岩应力的影响是一个比较长的时间过程，所以分析某年寒冷季节11月到翌年3月共5个月的隧洞洞口部位与隧洞中间部位自然通风150 d的温度-应力耦合分布规律，取每3 d的平均温度为1个温度荷载，共50个温度荷载，作用在洞壁围岩和进出口横断面上。模型上部施加地应力荷载且无约束，前后面y方向约束，左右面x方向约束，下部铰链约束。模型网格采用C3D8T八结点热耦合六面体单元，考虑到横断

面有圆形隧洞，给隧洞进行边布种，使得模型网格划分更合理。模型网格划分见图1所示。根据工程基本资料得到围岩力学参数，见表1。

表1 围岩力学参数

图1 模型网格划分

1.3 初始条件与边界条件

本文中围岩初始温度9℃，隧洞为通风边界，空气和围岩发生对流换热，因此洞壁为对流换热边界条件，通风初始温度4℃。由文献[20]可知，对流换热与流体的状态、热力学性质及固体壁面的几何尺寸、粗糙度有关，并且利用理想化的半圆拱形模型试验拟合出对流换热系数公式。本文研究隧洞不同位置围岩应力受通风温度、风速的影响，为满足单一变量原则，模型中假定：①隧洞围岩是连续均质、各向同性的，自然风流为连续介质，风速大小沿程不变。②洞壁围岩的换热条件是一致的，沿隧洞轴向和径向不发生变化。

依据文献[20]拟合出的计算公式求对流换热系数h。风速为1 m/s时，h=39.96 W/(m2·℃)；风速为2 m/s时，h=59.61 W/(m2·℃)；风速为3 m/s 时，h=88.93 W/(m2·℃)；风速为4 m/s时，h=132.67 W/(m2·℃)。隧洞埋深均取300 m，因为隧洞不同位置埋深相差不大，取最大极限有利于安全稳定，还能避免埋深变化的影响。

1.4 温度-应力耦合机理及求解方法

由文献[21]可知，温度-应力耦合有2个方面：一方面，是温度场对应力场的影响，围岩温度变化由温差产生的温度应力对围岩应力场的影响以及温度变化对围岩热力学参数的影响；另一方面，是应力场对温度场的影响，围岩应力的变化使围岩的内部结构发生变化，影响围岩的导热系数，进而影响围岩的温度场，与此同时，围岩内部变形会产生一定热量也会影响围岩温度场。

力学耦合机理：温度场和应力场之间通过某种力学作用相互耦合；参数耦合机理：温度场对岩体的物理性质参数的影响及岩体在不同风温、不同风速下热物理性质的改变影响温度场而表现出耦合作用。由文献[22]可知，温度场的变化对应力场的影响较大，而应力场在弹性变形范围内对温度场的影响非常小，从而可以忽略应力场对温度场的影响。本文为简化耦合计算，只考虑温度场对应力场的影响，忽略应力场改变对温度场的影响。

2 基于对流-导热耦合模型的围岩温度-应力分析

2.1 隧洞洞口和洞中位置温度模拟结果

图2、3分别为基于对流-导热耦合模型通风风速v为2 m/s时，距离隧洞洞口1 m和隧洞洞中位置寒冷季节不同天数的温度-应力耦合温度云图。表2为风速2 m/s时，自然通风情况下距离洞口1 m和隧洞洞中位置洞壁围岩的温度值。

表2 不同通风天数不同位置洞壁围岩温度 ℃

从图2、3可知，自然通风时间增加，对流-导热对隧洞周围围岩温度的影响范围逐渐增大，最终稳定在以隧洞中心为圆心半径10 m左右的范围，距离洞壁围岩2 m范围内的围岩温度变化最大，距离洞壁围岩越远，围岩的温度变化越小。从表2可知，隧洞60 d与30 d洞口的温差为1.15 ℃，洞中为1.31 ℃；隧洞90 d与30 d洞口的温差为6.26 ℃，洞中为6.48 ℃；隧洞120 d与30 d洞口的温差为1.88 ℃，洞中为2.48 ℃；隧洞150 d与30 d洞口的温差为9.01 ℃，洞中为9.425 ℃。隧洞在自然通风时，隧洞洞中位置产生的温差大于洞口。

图2 不同通风天数洞口位置围岩温度-应力耦合下温度分布(单位：℃)

2.2 隧洞洞口位置应力及塑性模拟结果

图4、5为通风风速2 m/s时，距离隧洞洞口1 m 的位置寒冷季节不同天数的温度-应力耦合最大主应力和最小主应力分布。表3为风速2 m/s时，距离洞口1 m位置围岩不同部位的主应力(拉应力为正，压应力为负)。

从图4、5和表3可知，顶拱、腰拱、底拱部位随着时间的增加主应力逐渐减小，主应力经过削减期到达稳定期。主要原因是围岩发生塑性变形使围岩内的能量释放，次要原因是低温下由温度差产生的温度应力(围岩单元受低温影响体积收缩产生拉应力)抵消了一部分围岩压应力。最大主应力与最小主应力差值在腰拱部位最大，并且随着时间的增加差值逐渐减小，最大差值为15.05 MPa。

图3 不同通风天数洞中位置围岩温度-应力耦合下温度分布(单位：℃)

图4 不同通风天数洞口位置温度-应力耦合下最大主应力分布(单位：Pa)

表3 不同天数围岩不同部位的主应力 MPa

提取图5中的数据，作出腰拱部位不同时间不同深度的最小主应力分布，见图6。从图6可知，腰拱部位随着时间的增加最小主应力逐渐减小，随着围岩深度的增加最小主应力逐渐减小至稳定，表明隧洞围岩经过应力削弱达到应力稳定。隧洞壁面围岩的最小主应力最大，这是由于洞壁处无衬砌结构并且出现应力集中现象造成的。

图5 不同通风天数洞口位置温度-应力耦合下最小主应力分布(单位：Pa)

图6 腰拱部位在不同天数不同深度的最小主应力变化

由文献[23]可知，在非均匀应力场测压力系数小于1时，水平方向的塑性区半径大于垂直方向的塑性区，塑性区形状近似为横蝶形。本文侧压力系数通过泊松比计算为0.43。图7为通风风速2 m/s时，距离隧洞洞口1 m的位置寒冷季节不同天数温度-应力耦合塑性应变值云图。从图7可知，塑性区形状大致为横蝶形，塑性应变最先出现在腰拱位置，并且随着通风时间的增加塑性应变值逐渐增大，塑性区逐渐由腰拱向上下两端延伸，塑性区范围逐渐增大。通风30 d时，塑性应变值最大出现在腰拱位置，为7.906×10-4；通风150 d，围岩塑性应变值最大出现在腰拱，为0.019 14。通风30 d时，塑性区半径约为1倍洞径，厚度约为0.69 m。随着通风时间的增加，塑性区向上下两端延伸，150 d时塑性区半径有微小变化，厚度变为2.08 m。

图7 不同通风天数洞口位置温度-应力耦合围岩塑性应变云图

2.3 隧洞洞中位置应力及塑性模拟结果

图8、9为通风风速2 m/s时隧洞洞中部位寒冷季节不同天数的温度-应力耦合最大主应力和最小主应力分布。通风150 d的最大主应力和最小主应力见表4。从图8、9可知，隧洞洞中部位的最大主应力与最小主应力变化规律和隧洞洞口位置基本一致。表4中应力值有明显的变化，洞口位置与洞中位置压应力值最大相差0.01 MPa。图2、3温度场结果可以很好解释这种变化规律，隧洞洞中位置温度高于洞口，自然通风时洞口位置的温差小于洞中位置的温差，所以洞口位置围岩体积收缩产生的拉应力小于洞中位置，低温产生的温度应力抵消了一部分围岩的压应力，使洞中位置主应力的压应力小于洞口位置。

图8 不同通风天数洞中位置温度-应力耦合下最大主应力分布(单位：Pa)

图9 不同通风天数洞中位置温度-应力耦合下最小主应力分布(单位：Pa)

表4 通风150 d不同位置的主应力 MPa

图10为通风风速2 m/s时隧洞洞中位置寒冷季节不同天数的温度-应力耦合塑性应变云图。不同通风时间隧洞洞口与隧洞洞中的塑性应变值见表5。从图10可知，通风150 d，隧洞洞中与隧洞洞口的塑性应变云图变化规律基本一致，塑性应变都是30 d时最先出现在腰拱位置，并且随着通风时间的增加塑性应变值逐渐增大，塑性区逐渐由腰拱向上下两端延伸，厚度逐渐增加。从表5可知，隧洞洞口位置的塑性应变大于隧洞洞中位置，一是，洞中位置的温差大，产生的温度应力大于洞口位置；二是，洞口和洞中位置的埋深一样，洞中位置的压应力小于洞口，故洞中位置的塑性应变也小。

图10 不同通风天数洞中位置温度-应力耦合围岩塑性应变云图

表5 隧洞洞口与洞中位置塑性应变 10-3

2.4 隧洞洞中位置不同风速下应力及塑性模拟结果

图11为隧洞洞中部位不同通风风速下经过寒冷季节150 d的温度-应力耦合最小主应力分布。表6为不同情况下温度-应力耦合的最小主应力值。从图11和表6可知，对照组(v=1 m/s且无温度荷载)与v=1 m/s的最大压应力差为0.16 MPa，最小压应力差为0.14 MPa，由此可以得出，通风150 d温度荷载在腰拱位置产生的最大温度应力(拉应力)为0.16 MPa，并且抵消了部分围岩压应力。随着通风风速的增加，最小压应力在缓慢减小，最大压应力也有缓慢减小，但表6中没体现(云图保留3位小数)。腰拱位置的压应力比温度应力大得多，故腰拱位置的应力受风速影响比顶拱和底拱受风速的影响小，文献[2]温度场结果可以很好地解释这种变化规律。当风速增加时，对流换热系数增加，隧洞洞壁温度逐渐降低，空气对流产生的温差逐渐增大，产生的温度应力(拉应力)逐渐增加从而抵消了一部分围岩压应力，因而风速增加围岩压应力在减小。

表6 隧洞洞中不同情况下最小主应力 MPa

图11 隧洞洞中不同风速作用下温度-应力耦合最小主应力分布(单位：Pa)

图12为隧洞洞中部位不同通风风速下经过寒冷季节150 d的温度-应力耦合塑性应变云图。表7为不同风速作用下温度-应力耦合的最大塑性应变值。从图12和表7可知，对照组(无温度荷载)的塑性应变值最大，为0.019 15，与v=1 m/s的塑性应变值相差0.000 41；随着通风风速的增加，对流换热系数增加，温度应力逐渐增大，抵消了围岩压应力，使围岩的塑性应变逐渐减小。由此可以得出，温度荷载对隧洞围岩的塑性应变影响较大，通风风速增大也可以影响隧洞围岩的塑性应变值，只是影响较小。