考慮熱交換的引水隧洞施工通風(fēng)兩相流模擬

王曉玲，劉 震，楊安林，張愛麗

(1. 天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072；2. 雅礱江流域水電開發(fā)有限公司，成都 610000)

考慮熱交換的引水隧洞施工通風(fēng)兩相流模擬

王曉玲1，劉 震1，楊安林2，張愛麗1

(1. 天津大學(xué)水利工程仿真與安全國家重點實驗室，天津 300072；2. 雅礱江流域水電開發(fā)有限公司，成都 610000)

針對深埋引水隧洞的儲熱作用所導(dǎo)致的施工通風(fēng)散煙困難問題，建立綜合考慮氣固熱交換、塵粒間相互耦合作用的三維非穩(wěn)態(tài)歐拉-拉格朗日兩相流紊流模型．塵粒間相互耦合作用力包括相間曳力、浮力、熱泳力、升力和虛擬質(zhì)量力等．結(jié)合某大型水電站深埋長距離引水隧洞施工通風(fēng)，探討隧洞內(nèi)含塵氣流的對流熱交換規(guī)律，研究隧洞內(nèi)風(fēng)流微觀運動機理以及粉塵運移規(guī)律．模擬結(jié)果表明：風(fēng)流在風(fēng)管出口下側(cè)形成渦流區(qū)，1號和2號引水隧洞風(fēng)流分布基本一致；工作面附近通風(fēng)降溫效果明顯，隧洞壁面的表面平均換熱系數(shù)沿工作面到風(fēng)管出口處呈先增大后減小的趨勢，表面平均換熱系數(shù)模擬值與經(jīng)驗公式計算結(jié)果基本吻合；施工通風(fēng)1,350,s左右時，粉塵濃度滿足規(guī)范要求．

引水隧洞；施工通風(fēng)兩相流；氣固熱交換；風(fēng)流運動；粉塵運移

大型深埋長距離引水隧洞工程一般具有高埋深、大洞徑、洞線長的特點．在這些大型深埋長引水隧洞工程的開挖中，廣泛采用鉆爆開挖這種經(jīng)濟高效的手段．由于鉆爆法產(chǎn)生大量的熱量，同時深埋長隧洞又具有儲熱作用，因此導(dǎo)致深埋引水隧洞施工通風(fēng)散煙困難．施工通風(fēng)排煙是確保施工進度順利完成的重要因素，它不僅關(guān)系到改善大型深埋長引水隧洞的施工環(huán)境，更是關(guān)系到施工人員的生命安全．因此，大型深埋長引水隧洞的施工通風(fēng)問題是一項具有挑戰(zhàn)性、亟待研究解決的關(guān)鍵科學(xué)問題．

目前國內(nèi)外對地下工程施工通風(fēng)的兩相流數(shù)值模擬主要集中于對礦井巷道的通風(fēng)設(shè)計方面. Laslandes等[1]和Klemens等[2]分別建立了二維Eulerian兩相流模型，研究了礦井巷道內(nèi)障礙物對礦塵分布的影響．楊勝來等[3]采用滑移擴散模型，模擬了礦井綜采工作面內(nèi)的風(fēng)速和粉塵濃度分布．陳贊成等[4]建立了高寒礦井穿脈巷道掘進的網(wǎng)格模型，利用Fluent軟件對不同掘進進尺和風(fēng)筒口距掘進工作面的距離在不同工況條件下的炮煙擴散規(guī)律進行了數(shù)值模擬，而在引水隧洞方面的兩相流模擬少有涉及．張靜等[5]采用了三維非穩(wěn)態(tài)歐拉兩相流模型模擬了云南南汀河引水隧洞獨頭掘進工作面爆破通風(fēng)過程，分析了隧洞內(nèi)的風(fēng)流結(jié)構(gòu)，研究了渦流、壁面粗糙度對粉塵擴散的影響．

綜上所述，在以上地下工程施工通風(fēng)兩相流模擬研究中均沒有考慮氣固熱交換對施工通風(fēng)的影響.為此，筆者建立了考慮熱交換影響的引水隧洞施工通風(fēng)三維非穩(wěn)態(tài)歐拉-拉格朗日兩相流紊流模型，在動量方程中考慮浮力、相間曳力、熱泳力、升力和虛擬質(zhì)量力的作用，結(jié)合某大型水電站引水隧洞施工通風(fēng)進行了數(shù)值模擬，分析了隧洞群內(nèi)風(fēng)流場和溫度場分布、熱交換以及粉塵運移規(guī)律，以期為工程應(yīng)用提供理論依據(jù)和參考．

1 數(shù)學(xué)模型

采用歐拉-拉格朗日模型，氣相采用標(biāo)準(zhǔn)k-ε 紊流模型，顆粒相采用離散相模型(DPM)，粉塵顆粒運動軌跡的追蹤采用考慮顆粒湍流擴散影響的隨機顆粒軌道模型．動量方程考慮溫度差和密度差引起的浮力、相間曳力、熱泳力、升力和虛擬質(zhì)量力的作用. 采用有限體積法對控制方程進行離散，采用PISO算法進行求解[6]．

1.1 控制方程

1) 氣相控制方程

選取標(biāo)準(zhǔn)k-ε紊流模型為氣相的紊流封閉方程．采用Boussinesq近似，除動量方程的浮力計算外，其余方程的氣體密度為常數(shù)．

氣相連續(xù)性方程為

式中：t為時間，s；cρ為氣體密度，3kg/m；cu為氣體速度矢量，m/s．

氣相動量方程為

標(biāo)準(zhǔn)k-ε紊流模型為

式中：ui為時均速度；tμ為紊動黏性系數(shù)μ為分子動力黏性系數(shù)，(N·m)/s；ck為紊動動能，m2/s2；cε為紊動能耗散率，m2/s3；kG為由平均速度梯度產(chǎn)生的紊動動能；kS、Sε分別為k動量方程和ε輸運方程的源項；bG為由于浮力引起的紊動動能k的產(chǎn)生項；1Cε、2Cε、3Cε、Cμ、kσ和εσ為經(jīng)驗常數(shù)，C1ε=1.44，C2ε=1.92，C3ε=0，Cμ=0.09，σk=1.0，σε=1.3．

2) 能量守恒方程

以溫度T為變量的能量守恒方程[7]為

式中：pc為比定壓熱容；T為溫度；K為流體的傳熱系數(shù)；TS為黏性耗散項，包括流體的內(nèi)熱源和由于黏性作用流體機械能轉(zhuǎn)換為熱能的部分．

3) 顆粒相控制方程

考慮湍流脈動速度對離散相顆粒軌跡的影響，采用隨機軌道模型計算顆粒與湍流的相互作用．顆粒離散相的控制方程為

式中：su為顆粒相的瞬時速度，m/s；sm為顆粒相質(zhì)量，kg；DS為氣相和顆粒相間由密度差、溫度差引起的浮力；sS為氣相對顆粒相的作用力，包括曳力DF、熱泳力TF、壓強梯度力PF、虛擬質(zhì)量力VMF和升力FL；Fss為顆粒間碰撞力；CD為曳力系數(shù)，取0.44；ds為顆粒直徑，m；DT,s為熱泳力系數(shù)；Vs為顆粒體積，m3；?p為氣相的壓強梯度，N/m3；νs′為顆粒相的脈動速度，m/s．

1.2 邊界條件

假定流體為不可壓縮、非穩(wěn)態(tài)紊流，滿足Boussinesq假設(shè)，風(fēng)管進口風(fēng)速分布均勻，模擬邊界條件如下．

1) 進口邊界

uin=Q/S，S為風(fēng)管進口斷面面積，m2；Q為通風(fēng)風(fēng)量，m3/s；進口速度分量in0v=，in0w=；湍流動能k和紊流動能耗散率ε由經(jīng)驗公式計算得到[8]．

2) 出口邊界

采用壓力出口邊界條件，隧洞出口與外界大氣相通，出口處表壓設(shè)為0，溫度設(shè)為外界溫度．

3) 固體壁面

隧洞面和工作面均為固體壁面，按固壁定律處理，采用壁面無滑移條件，流體計算域邊界采用考慮粗糙度影響的壁面函數(shù)法[9]．

4) 壁面熱邊界條件

采用對流換熱邊界條件，且引入溫度壁面函數(shù)法計算求解換熱系數(shù)，即

5) 顆粒相邊界條件

粉塵顆粒源設(shè)在工作面，除隧洞出口邊界設(shè)為逃逸外，其余的壁面均設(shè)為反射，并且假設(shè)顆粒碰撞為完全彈性碰撞．

1.3 初始條件

考慮某水利水電地下工程引水隧洞的圍巖主要為碳酸鹽巖及少量砂巖、板巖、綠泥石片巖等，故取干放炮的粉塵產(chǎn)量為400,mg/m3，假設(shè)進口處的粉塵濃度為0．

圖1 引水隧洞掘進至550,m時物理模型(單位：m)Fig.1 Diversion tunnel model with 550,m tunneling length(unit：m)

2 引水隧洞施工通風(fēng)兩相流數(shù)值模擬分析

某大型水電站工程樞紐主要由首部攔河閘、引水系統(tǒng)和尾部地下廠房3部分組成，其引水隧洞洞群穿越主體山脈，上覆巖體一般埋深1,500～2,000,m，最大埋深約為2,525,m，具有埋深大、洞線長、洞徑大的特點．

以引水隧洞群1號、2號引水隧洞掘進至550,m時物理模型為例，研究引水隧洞內(nèi)風(fēng)流場、溫度場分布以及粉塵運移規(guī)律．引水隧洞開挖斷面直徑為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(寬×高)的城門洞形；風(fēng)管直徑2,m，置于隧洞頂部，風(fēng)管出風(fēng)口距工作面(A為隧洞斷面面積)．采用非結(jié)構(gòu)化貼體網(wǎng)格劃分技術(shù)和局部網(wǎng)格加密技術(shù)建立網(wǎng)格模型，主體區(qū)域網(wǎng)格尺寸為1,m×1,m×1,m，網(wǎng)格總數(shù)為465,633個，其物理模型和計算網(wǎng)格如圖1和圖2所示．

圖2 引水隧洞掘進至550,m時計算網(wǎng)格Fig.2 Diversion tunnel mesh model with 550,m tunneling length

2.1 風(fēng)流場模擬結(jié)果分析

圖3(a)和(b)分別為通風(fēng)50,s時1號工作面和2號工作面附近風(fēng)流場分布剖面．由圖3可知，通風(fēng)50,s時，風(fēng)流從風(fēng)管射出沿隧洞頂壁到達工作面，速度逐漸減小，在隧洞頂部形成射流區(qū)，見圖3(a)Ⅱ區(qū)和圖3(b)Ⅱ區(qū)．由于風(fēng)流的連續(xù)性和掘進空間的局限性，風(fēng)流到達工作面后沿工作面形成貼壁射流區(qū)，見圖3(a)Ⅲ區(qū)和圖3(b)Ⅲ區(qū)；并沿隧洞底部形成回流區(qū)，見圖3(a)Ⅳ區(qū)和圖3(b)Ⅳ區(qū)．在射流區(qū)和回流區(qū)之間，風(fēng)流形成了渦流分布，見圖3(a)Ⅰ區(qū)和圖3(b)Ⅰ區(qū)，但速度較小約為2,m/s．對比圖3(a)和(b)可知，兩引水隧洞風(fēng)流場分布趨勢基本一致．

2.2 熱交換模擬結(jié)果分析

圖4(a)和(b)分別為通風(fēng)初始時刻和通風(fēng)1,000,s時隧洞壁面溫度場分布，由圖4(a)可知，在通風(fēng)初始時刻，由于工作面處鉆爆施工產(chǎn)生大量熱量導(dǎo)致工作面處溫度升高，最高達到28,℃．由圖4(b)可知，隨著通風(fēng)時間的增加，風(fēng)流與隧洞壁面間的氣固熱交換不斷進行，工作面附近的壁面溫度逐漸降低，并維持在21,℃左右．

圖3 t=50,s時工作面附近風(fēng)流場分布剖面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

對流換熱系數(shù)是計算巷道壁面與風(fēng)流之間熱交換量的重要參數(shù)．本文采用溫度壁面函數(shù)法求解換熱系數(shù)，為驗證該方法可靠性，采用根據(jù)工程實測資料[10]所得出的計算公式計算的換熱系數(shù)進行對比分析．計算公式為

式中：λ為空氣導(dǎo)熱系數(shù)，W/(m2·K)；ed為斷面水力直徑，m；Re為雷諾數(shù)．

圖5為隧洞壁面熱交換系數(shù)分布．由圖5可知，工作面和風(fēng)管出口附近熱交換系數(shù)較大，這是由于工作面處的溫度比其他壁面處的溫度高，風(fēng)管出口附近風(fēng)流速度較其他位置處風(fēng)速大，故工作面和風(fēng)管出口附近熱交換強度均較大所致．圖6為隧洞壁面表面平均換熱系數(shù)模擬值和經(jīng)驗公式計算值對比．由圖6可知，隧洞壁面的表面平均換熱系數(shù)分布規(guī)律和經(jīng)驗公式計算值所得的分布規(guī)律基本一致，均呈現(xiàn)沿工作面至風(fēng)管出口(距工作面30,m處)呈先增大后減小的趨勢．但由于采用的經(jīng)驗公式是根據(jù)某個地下實際工程得出，故與本工程存在一定差異，其中平均誤差為7.8%，最大誤差為14.3%．

圖5 隧洞壁面熱交換系數(shù)分布Fig.5 Distribution of heat transfer coefficient on the wall

圖6 表面平均換熱系數(shù)對比Fig.6 Comparison of surface average heat transfer coefficient

2.3 粉塵運移模擬結(jié)果分析

圖7為風(fēng)流場穩(wěn)定時的粉塵顆粒軌跡及停留時間分布．由圖7可知，粉塵在1號工作面和2號工作面附近停留時間較短，為100,s左右，在風(fēng)流的攜帶作用下，粉塵由工作面不斷向支洞方向運移并最終從出口排出，粉塵顆粒在隧洞內(nèi)的最大停留時間為2,000,s左右．由圖中局部放大處可知，粉塵基本不在1號橫通道和支洞連接處停留，說明1號和2號引水隧洞間基本沒有粉塵顆粒通過，1號工作面處產(chǎn)生的粉塵主要通過引1號支洞由出口1排出，2號工作面處產(chǎn)生的粉塵主要通過引2號支洞由出口2排出．

引水隧洞施工通風(fēng)的排塵率累積曲線見圖8．施工通風(fēng)約960,s時，排塵率達86.1%左右，即86.1%的粉塵從隧洞出口處排出；施工通風(fēng)1,350,s時，排塵率為97.7%左右，隧洞出口的粉塵質(zhì)量流率為4.12× 10-5,kg/s，粉塵質(zhì)量濃度為0.74,mg/m3，小于規(guī)范允許值2.00,mg/m3，滿足隧洞施工通風(fēng)環(huán)境要求．

圖7 粉塵顆粒軌跡及停留時間分布Fig.7 Dust particle tracks and distribution of dust residence time

圖8 排塵率累積曲線Fig.8 Accumulation curve of dust discharge rate

3 結(jié) 語

在考慮氣固熱交換影響的情況下，建立了三維非穩(wěn)態(tài)歐拉-拉格朗日兩相流紊流模型，并模擬了某水電站深埋引水隧洞施工通風(fēng)過程，分析了隧洞內(nèi)風(fēng)流場和溫度場分布以及粉塵運移規(guī)律，探討了傳熱分析的關(guān)鍵參數(shù)換熱系數(shù)的分布規(guī)律，為深埋引水隧洞施工通風(fēng)提供熱工計算依據(jù)．

模擬結(jié)果表明：風(fēng)流從風(fēng)管射出沿隧洞頂部形成射流區(qū)，到達工作面后形成貼壁射流區(qū)，并沿隧洞底部形成回流區(qū)；射流區(qū)和回流區(qū)之間，風(fēng)流形成渦流分布，且1號和2號引水隧洞風(fēng)流場分布趨勢基本一致；隨著通風(fēng)時間的增加，工作面附近溫度逐漸降低，基本達到了通風(fēng)降溫效果；隧洞壁面的表面平均換熱系數(shù)沿工作面到風(fēng)管出口處呈先增大后減小的趨勢，隧洞壁面的表面平均換熱系數(shù)模擬值與經(jīng)驗公式計算結(jié)果基本吻合；施工通風(fēng)1,350,s左右時，粉塵濃度滿足規(guī)范要求．

［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］ 楊勝來，黃元平. 綜采工作面粉塵濃度分布的數(shù)值解法[J]. 中國安全科學(xué)學(xué)報，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］ 陳贊成，楊 鵬，呂文生，等. 高寒礦井穿脈巷道掘進炮煙擴散規(guī)律的數(shù)值模擬[J]. 北京科技大學(xué)學(xué)報，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］ 張 靜，王曉玲，陳紅超，等. 引水隧洞獨頭掘進工作面風(fēng)流和粉塵擴散的模擬[J]. 水力發(fā)電學(xué)報，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］ 王福軍. 計算機流體動力學(xué)分析[M]. 北京：清華大學(xué)出版社，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］ 牟靈泉. 地道風(fēng)降溫計算與應(yīng)用[M]. 北京：中國建筑工業(yè)出版社，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.

國家自然科學(xué)基金資助項目(51179121，51209160)；國家自然科學(xué)基金創(chuàng)新研究群體科學(xué)基金資助項目(51021004)．

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

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