考虑热交换的引水隧洞施工通风两相流模拟

王晓玲，刘 震，杨安林，张爱丽

(1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 雅砻江流域水电开发有限公司，成都 610000)

考虑热交换的引水隧洞施工通风两相流模拟

王晓玲1，刘 震1，杨安林2，张爱丽1

(1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 雅砻江流域水电开发有限公司，成都 610000)

针对深埋引水隧洞的储热作用所导致的施工通风散烟困难问题，建立综合考虑气固热交换、尘粒间相互耦合作用的三维非稳态欧拉-拉格朗日两相流紊流模型．尘粒间相互耦合作用力包括相间曳力、浮力、热泳力、升力和虚拟质量力等．结合某大型水电站深埋长距离引水隧洞施工通风，探讨隧洞内含尘气流的对流热交换规律，研究隧洞内风流微观运动机理以及粉尘运移规律．模拟结果表明：风流在风管出口下侧形成涡流区，1号和2号引水隧洞风流分布基本一致；工作面附近通风降温效果明显，隧洞壁面的表面平均换热系数沿工作面到风管出口处呈先增大后减小的趋势，表面平均换热系数模拟值与经验公式计算结果基本吻合；施工通风1,350,s左右时，粉尘浓度满足规范要求．

引水隧洞；施工通风两相流；气固热交换；风流运动；粉尘运移

大型深埋长距离引水隧洞工程一般具有高埋深、大洞径、洞线长的特点．在这些大型深埋长引水隧洞工程的开挖中，广泛采用钻爆开挖这种经济高效的手段．由于钻爆法产生大量的热量，同时深埋长隧洞又具有储热作用，因此导致深埋引水隧洞施工通风散烟困难．施工通风排烟是确保施工进度顺利完成的重要因素，它不仅关系到改善大型深埋长引水隧洞的施工环境，更是关系到施工人员的生命安全．因此，大型深埋长引水隧洞的施工通风问题是一项具有挑战性、亟待研究解决的关键科学问题．

目前国内外对地下工程施工通风的两相流数值模拟主要集中于对矿井巷道的通风设计方面. Laslandes等[1]和Klemens等[2]分别建立了二维Eulerian两相流模型，研究了矿井巷道内障碍物对矿尘分布的影响．杨胜来等[3]采用滑移扩散模型，模拟了矿井综采工作面内的风速和粉尘浓度分布．陈赞成等[4]建立了高寒矿井穿脉巷道掘进的网格模型，利用Fluent软件对不同掘进进尺和风筒口距掘进工作面的距离在不同工况条件下的炮烟扩散规律进行了数值模拟，而在引水隧洞方面的两相流模拟少有涉及．张静等[5]采用了三维非稳态欧拉两相流模型模拟了云南南汀河引水隧洞独头掘进工作面爆破通风过程，分析了隧洞内的风流结构，研究了涡流、壁面粗糙度对粉尘扩散的影响．

综上所述，在以上地下工程施工通风两相流模拟研究中均没有考虑气固热交换对施工通风的影响.为此，笔者建立了考虑热交换影响的引水隧洞施工通风三维非稳态欧拉-拉格朗日两相流紊流模型，在动量方程中考虑浮力、相间曳力、热泳力、升力和虚拟质量力的作用，结合某大型水电站引水隧洞施工通风进行了数值模拟，分析了隧洞群内风流场和温度场分布、热交换以及粉尘运移规律，以期为工程应用提供理论依据和参考．

1 数学模型

采用欧拉-拉格朗日模型，气相采用标准k-ε 紊流模型，颗粒相采用离散相模型(DPM)，粉尘颗粒运动轨迹的追踪采用考虑颗粒湍流扩散影响的随机颗粒轨道模型．动量方程考虑温度差和密度差引起的浮力、相间曳力、热泳力、升力和虚拟质量力的作用. 采用有限体积法对控制方程进行离散，采用PISO算法进行求解[6]．

1.1 控制方程

1) 气相控制方程

选取标准k-ε紊流模型为气相的紊流封闭方程．采用Boussinesq近似，除动量方程的浮力计算外，其余方程的气体密度为常数．

气相连续性方程为

式中：t为时间，s；cρ为气体密度，3kg/m；cu为气体速度矢量，m/s．

气相动量方程为

标准k-ε紊流模型为

式中：ui为时均速度；tμ为紊动黏性系数μ为分子动力黏性系数，(N·m)/s；ck为紊动动能，m2/s2；cε为紊动能耗散率，m2/s3；kG为由平均速度梯度产生的紊动动能；kS、Sε分别为k动量方程和ε输运方程的源项；bG为由于浮力引起的紊动动能k的产生项；1Cε、2Cε、3Cε、Cμ、kσ和εσ为经验常数，C1ε=1.44，C2ε=1.92，C3ε=0，Cμ=0.09，σk=1.0，σε=1.3．

2) 能量守恒方程

以温度T为变量的能量守恒方程[7]为

式中：pc为比定压热容；T为温度；K为流体的传热系数；TS为黏性耗散项，包括流体的内热源和由于黏性作用流体机械能转换为热能的部分．

3) 颗粒相控制方程

考虑湍流脉动速度对离散相颗粒轨迹的影响，采用随机轨道模型计算颗粒与湍流的相互作用．颗粒离散相的控制方程为

式中：su为颗粒相的瞬时速度，m/s；sm为颗粒相质量，kg；DS为气相和颗粒相间由密度差、温度差引起的浮力；sS为气相对颗粒相的作用力，包括曳力DF、热泳力TF、压强梯度力PF、虚拟质量力VMF和升力FL；Fss为颗粒间碰撞力；CD为曳力系数，取0.44；ds为颗粒直径，m；DT,s为热泳力系数；Vs为颗粒体积，m3；∇p为气相的压强梯度，N/m3；νs′为颗粒相的脉动速度，m/s．

1.2 边界条件

假定流体为不可压缩、非稳态紊流，满足Boussinesq假设，风管进口风速分布均匀，模拟边界条件如下．

1) 进口边界

uin=Q/S，S为风管进口断面面积，m2；Q为通风风量，m3/s；进口速度分量in0v=，in0w=；湍流动能k和紊流动能耗散率ε由经验公式计算得到[8]．

2) 出口边界

采用压力出口边界条件，隧洞出口与外界大气相通，出口处表压设为0，温度设为外界温度．

3) 固体壁面

隧洞面和工作面均为固体壁面，按固壁定律处理，采用壁面无滑移条件，流体计算域边界采用考虑粗糙度影响的壁面函数法[9]．

4) 壁面热边界条件

采用对流换热边界条件，且引入温度壁面函数法计算求解换热系数，即

5) 颗粒相边界条件

粉尘颗粒源设在工作面，除隧洞出口边界设为逃逸外，其余的壁面均设为反射，并且假设颗粒碰撞为完全弹性碰撞．

1.3 初始条件

考虑某水利水电地下工程引水隧洞的围岩主要为碳酸盐岩及少量砂岩、板岩、绿泥石片岩等，故取干放炮的粉尘产量为400,mg/m3，假设进口处的粉尘浓度为0．

图1 引水隧洞掘进至550,m时物理模型(单位：m)Fig.1 Diversion tunnel model with 550,m tunneling length(unit：m)

2 引水隧洞施工通风两相流数值模拟分析

某大型水电站工程枢纽主要由首部拦河闸、引水系统和尾部地下厂房3部分组成，其引水隧洞洞群穿越主体山脉，上覆岩体一般埋深1,500～2,000,m，最大埋深约为2,525,m，具有埋深大、洞线长、洞径大的特点．

以引水隧洞群1号、2号引水隧洞掘进至550,m时物理模型为例，研究引水隧洞内风流场、温度场分布以及粉尘运移规律．引水隧洞开挖断面直径为13,m，断面为马蹄形，分上下2个断面开挖，上半部分先行开挖，开挖断面面积约为98,m2；引1号支洞长664,m，纵坡6.88%，断面为8.50,m×7.24,m(宽×高)的城门洞形；引2号支洞长820,m，纵坡6.96%，断面为8.50,m×7.24,m(宽×高)的城门洞形；1号横通道长47,m，断面为8.50,m×7.24,m(宽×高)的城门洞形；风管直径2,m，置于隧洞顶部，风管出风口距工作面(A为隧洞断面面积)．采用非结构化贴体网格划分技术和局部网格加密技术建立网格模型，主体区域网格尺寸为1,m×1,m×1,m，网格总数为465,633个，其物理模型和计算网格如图1和图2所示．

图2 引水隧洞掘进至550,m时计算网格Fig.2 Diversion tunnel mesh model with 550,m tunneling length

2.1 风流场模拟结果分析

图3(a)和(b)分别为通风50,s时1号工作面和2号工作面附近风流场分布剖面．由图3可知，通风50,s时，风流从风管射出沿隧洞顶壁到达工作面，速度逐渐减小，在隧洞顶部形成射流区，见图3(a)Ⅱ区和图3(b)Ⅱ区．由于风流的连续性和掘进空间的局限性，风流到达工作面后沿工作面形成贴壁射流区，见图3(a)Ⅲ区和图3(b)Ⅲ区；并沿隧洞底部形成回流区，见图3(a)Ⅳ区和图3(b)Ⅳ区．在射流区和回流区之间，风流形成了涡流分布，见图3(a)Ⅰ区和图3(b)Ⅰ区，但速度较小约为2,m/s．对比图3(a)和(b)可知，两引水隧洞风流场分布趋势基本一致．

2.2 热交换模拟结果分析

图4(a)和(b)分别为通风初始时刻和通风1,000,s时隧洞壁面温度场分布，由图4(a)可知，在通风初始时刻，由于工作面处钻爆施工产生大量热量导致工作面处温度升高，最高达到28,℃．由图4(b)可知，随着通风时间的增加，风流与隧洞壁面间的气固热交换不断进行，工作面附近的壁面温度逐渐降低，并维持在21,℃左右．

图3 t=50,s时工作面附近风流场分布剖面Fig.3 Distribution of airflow field near the working face after 50,s

图4 不同时刻下隧洞壁面温度场分布Fig.4 Distribution of temperature field on the wall at different times

对流换热系数是计算巷道壁面与风流之间热交换量的重要参数．本文采用温度壁面函数法求解换热系数，为验证该方法可靠性，采用根据工程实测资料[10]所得出的计算公式计算的换热系数进行对比分析．计算公式为

式中：λ为空气导热系数，W/(m2·K)；ed为断面水力直径，m；Re为雷诺数．

图5为隧洞壁面热交换系数分布．由图5可知，工作面和风管出口附近热交换系数较大，这是由于工作面处的温度比其他壁面处的温度高，风管出口附近风流速度较其他位置处风速大，故工作面和风管出口附近热交换强度均较大所致．图6为隧洞壁面表面平均换热系数模拟值和经验公式计算值对比．由图6可知，隧洞壁面的表面平均换热系数分布规律和经验公式计算值所得的分布规律基本一致，均呈现沿工作面至风管出口(距工作面30,m处)呈先增大后减小的趋势．但由于采用的经验公式是根据某个地下实际工程得出，故与本工程存在一定差异，其中平均误差为7.8%，最大误差为14.3%．

图5 隧洞壁面热交换系数分布Fig.5 Distribution of heat transfer coefficient on the wall

图6 表面平均换热系数对比Fig.6 Comparison of surface average heat transfer coefficient

2.3 粉尘运移模拟结果分析

图7为风流场稳定时的粉尘颗粒轨迹及停留时间分布．由图7可知，粉尘在1号工作面和2号工作面附近停留时间较短，为100,s左右，在风流的携带作用下，粉尘由工作面不断向支洞方向运移并最终从出口排出，粉尘颗粒在隧洞内的最大停留时间为2,000,s左右．由图中局部放大处可知，粉尘基本不在1号横通道和支洞连接处停留，说明1号和2号引水隧洞间基本没有粉尘颗粒通过，1号工作面处产生的粉尘主要通过引1号支洞由出口1排出，2号工作面处产生的粉尘主要通过引2号支洞由出口2排出．

引水隧洞施工通风的排尘率累积曲线见图8．施工通风约960,s时，排尘率达86.1%左右，即86.1%的粉尘从隧洞出口处排出；施工通风1,350,s时，排尘率为97.7%左右，隧洞出口的粉尘质量流率为4.12× 10-5,kg/s，粉尘质量浓度为0.74,mg/m3，小于规范允许值2.00,mg/m3，满足隧洞施工通风环境要求．

图7 粉尘颗粒轨迹及停留时间分布Fig.7 Dust particle tracks and distribution of dust residence time

图8 排尘率累积曲线Fig.8 Accumulation curve of dust discharge rate

3 结 语

在考虑气固热交换影响的情况下，建立了三维非稳态欧拉-拉格朗日两相流紊流模型，并模拟了某水电站深埋引水隧洞施工通风过程，分析了隧洞内风流场和温度场分布以及粉尘运移规律，探讨了传热分析的关键参数换热系数的分布规律，为深埋引水隧洞施工通风提供热工计算依据．

模拟结果表明：风流从风管射出沿隧洞顶部形成射流区，到达工作面后形成贴壁射流区，并沿隧洞底部形成回流区；射流区和回流区之间，风流形成涡流分布，且1号和2号引水隧洞风流场分布趋势基本一致；随着通风时间的增加，工作面附近温度逐渐降低，基本达到了通风降温效果；隧洞壁面的表面平均换热系数沿工作面到风管出口处呈先增大后减小的趋势，隧洞壁面的表面平均换热系数模拟值与经验公式计算结果基本吻合；施工通风1,350,s左右时，粉尘浓度满足规范要求．

［1］ Laslandes S，Sacé C. Transport of particles by a turbulent flow around an obstacle—A numerical and a wind tunnel approach [J]. Journal of Wind Engineering and Industrial Aerodynamics，1998，74/75/76：577-587.

［2］ Klemens R，Kosinski P，Wolanski P，et al. Numerical study of dust lifting in a channel with vertical obstacles [J]. Journal of Loss Prevention in the Process Industries，2001，14(6)：469-473.

［3］ 杨胜来，黄元平. 综采工作面粉尘浓度分布的数值解法[J]. 中国安全科学学报，1996，6(增)：64-68.

Yang Shenglai，Huang Yuanping. Numerical simulation of dust distribution in longwall face in coal mine[J]. China Safety Science Journal，1996，6(Suppl)：64-68(in Chinese).

［4］ 陈赞成，杨 鹏，吕文生，等. 高寒矿井穿脉巷道掘进炮烟扩散规律的数值模拟[J]. 北京科技大学学报，2011，33(5)：521-525.

Chen Zancheng，Yang Peng，Lü Wensheng，et al. Numerical simulation on the diffusion law of blasting fume during roadway tunneling across a vein in an alpine mine[J]. Journal of University of Science and Technology Beijing，2011，33(5)：521-525(in Chinese).

［5］ 张 静，王晓玲，陈红超，等. 引水隧洞独头掘进工作面风流和粉尘扩散的模拟[J]. 水力发电学报，2008，27(1)：111-117.

Zhang Jing，Wang Xiaoling，Chen Hongchao，et al. Simulation on ventilation and dust diffusion on heading face of the diversion tunnel[J]. Journal of Hydroelectric Engineering，2008，27(1)：111-117(in Chinese).

［6］ 王福军. 计算机流体动力学分析[M]. 北京：清华大学出版社，2004.

Wang Fujun. Analysis of Computational Fluid Dynamics[M]. Beijing：Tsinghua University Press，2004(in Chinese).

［7］ Issa R I. Solution of implicitly discretized fluid flow equations by operator splitting[J]. Comput Phys，1986，62(1)：40-65.

［8］ Launder B E，Spalding D B. The numerical computation of turbulent flows[J]. Comput Method Appl Mech Eng，1974，3(2)：269-289.

［9］ Lacasse D，Turgeon E，Pelletier D. On the judicious use of the k-ε model，wall functions and adaptivity[J]. International Journal of Thermal Sciences，2004，43(10)：925-938.

［10］ 牟灵泉. 地道风降温计算与应用[M]. 北京：中国建筑工业出版社，1982.

Mou Lingquan. Computation and Application of Air Through Tunnel Cooling[M]. Beijing：China Architecture and Building Press，1982(in Chinese).

Simulation for Construction Ventilation Two-Phase Flow in Diversion Tunnel Considering Heat Exchange

Wang Xiaoling1，Liu Zhen1，Yang Anlin2，Zhang Aili1
(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2. Yalong River Hydropower Development Company，Chengdu 610000，China)

Smoke extraction difficulties caused by storage heat effects exist in the construction ventilation of deep diversion tunnel. A three-dimensional unsteady Euler-Lagrange two-phase flow turbulence model considering the gassolid heat exchange and the intercoupling between dust particles was proposed. The intercoupling between dust particles included the effects of buoyancy，inter-phase drag，thermophoretic force，lift force and virtual mass force. The construction ventilation of a large-scale hydropower long-distance diversion tunnel was used as a case. The dusty air convection heat transfer law in the tunnel was discussed. The microscopic airflow movement mechanism and the dust migration law were analyzed. Results show that an eddy zone is formed under the duct outlet，and that the airflow distribution tendency of No.1 diversion tunnel is similar to that of No.2 diversion tunnel. With the ventilation time increasing，the temperature near the heading face apparently drops. The surface averageheat transfer coefficients of the tunnel have a tendency to increase firstly and then diminish from the heading face to the duct outlet. The simulation values of heat transfer coefficients are in good agreement with the experiential values. Dust concentration meets the standard after about 1,350,s.

diversion tunnel；construction ventilation two-phase flow；gas-solid heat exchange；airflow movement；dust migration

TV554

A

0493-2137(2013)04-0322-06

DOI 10.11784/tdxb20130406

2012-11-26；

2013-01-04.

国家自然科学基金资助项目(51179121，51209160)；国家自然科学基金创新研究群体科学基金资助项目(51021004)．

王晓玲（1968— ），女，博士，教授.

王晓玲，wangxl @tju.edu.cn.