王曉玲，劉 震，楊安林，張愛(ài)麗

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

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

王曉玲1，劉 震1，楊安林2，張愛(ài)麗1

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

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

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

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

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

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

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

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

1.1 控制方程

1) 氣相控制方程

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

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

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

氣相動(dòng)量方程為

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

式中：ui為時(shí)均速度；tμ為紊動(dòng)黏性系數(shù)μ為分子動(dòng)力黏性系數(shù)，(N·m)/s；ck為紊動(dòng)動(dòng)能，m2/s2；cε為紊動(dòng)能耗散率，m2/s3；kG為由平均速度梯度產(chǎn)生的紊動(dòng)動(dòng)能；kS、Sε分別為k動(dòng)量方程和ε輸運(yùn)方程的源項(xiàng)；bG為由于浮力引起的紊動(dòng)動(dòng)能k的產(chǎn)生項(xiàng)；1Cε、2Cε、3Cε、Cμ、kσ和εσ為經(jīng)驗(yàn)常數(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為黏性耗散項(xiàng)，包括流體的內(nèi)熱源和由于黏性作用流體機(jī)械能轉(zhuǎn)換為熱能的部分．

3) 顆粒相控制方程

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

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

1.2 邊界條件

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

1) 進(jìn)口邊界

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

2) 出口邊界

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

3) 固體壁面

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

4) 壁面熱邊界條件

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

5) 顆粒相邊界條件

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

1.3 初始條件

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

圖1 引水隧洞掘進(jìn)至550,m時(shí)物理模型(單位：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，具有埋深大、洞線長(zhǎng)、洞徑大的特點(diǎn)．

以引水隧洞群1號(hào)、2號(hào)引水隧洞掘進(jìn)至550,m時(shí)物理模型為例，研究引水隧洞內(nèi)風(fēng)流場(chǎng)、溫度場(chǎng)分布以及粉塵運(yùn)移規(guī)律．引水隧洞開(kāi)挖斷面直徑為13,m，斷面為馬蹄形，分上下2個(gè)斷面開(kāi)挖，上半部分先行開(kāi)挖，開(kāi)挖斷面面積約為98,m2；引1號(hào)支洞長(zhǎng)664,m，縱坡6.88%，斷面為8.50,m×7.24,m(寬×高)的城門洞形；引2號(hào)支洞長(zhǎng)820,m，縱坡6.96%，斷面為8.50,m×7.24,m(寬×高)的城門洞形；1號(hào)橫通道長(zhǎng)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個(gè)，其物理模型和計(jì)算網(wǎng)格如圖1和圖2所示．

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

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

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

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

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

圖3 t=50,s時(shí)工作面附近風(fēng)流場(chǎng)分布剖面Fig.3 Distribution of airflow field near the working face after 50,s

圖4 不同時(shí)刻下隧洞壁面溫度場(chǎng)分布Fig.4 Distribution of temperature field on the wall at different times

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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

模擬結(jié)果表明：風(fēng)流從風(fēng)管射出沿隧洞頂部形成射流區(qū)，到達(dá)工作面后形成貼壁射流區(qū)，并沿隧洞底部形成回流區(qū)；射流區(qū)和回流區(qū)之間，風(fēng)流形成渦流分布，且1號(hào)和2號(hào)引水隧洞風(fēng)流場(chǎng)分布趨勢(shì)基本一致；隨著通風(fēng)時(shí)間的增加，工作面附近溫度逐漸降低，基本達(dá)到了通風(fēng)降溫效果；隧洞壁面的表面平均換熱系數(shù)沿工作面到風(fēng)管出口處呈先增大后減小的趨勢(shì)，隧洞壁面的表面平均換熱系數(shù)模擬值與經(jīng)驗(yàn)公式計(jì)算結(jié)果基本吻合；施工通風(fēng)1,350,s左右時(shí)，粉塵濃度滿足規(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］ 楊勝來(lái)，黃元平. 綜采工作面粉塵濃度分布的數(shù)值解法[J]. 中國(guó)安全科學(xué)學(xué)報(bào)，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ìn)炮煙擴(kuò)散規(guī)律的數(shù)值模擬[J]. 北京科技大學(xué)學(xué)報(bào)，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］ 張 靜，王曉玲，陳紅超，等. 引水隧洞獨(dú)頭掘進(jìn)工作面風(fēng)流和粉塵擴(kuò)散的模擬[J]. 水力發(fā)電學(xué)報(bào)，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］ 王福軍. 計(jì)算機(jī)流體動(dòng)力學(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)降溫計(jì)算與應(yīng)用[M]. 北京：中國(guó)建筑工業(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.

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

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

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