陳麗萍，程璟濤，蔣軍成，鄧廣發(fā)

(1.南京工業(yè)大學(xué)城市建設(shè)與安全工程學(xué)院，江蘇南京 210009;2.江蘇方天電力技術(shù)有限公司，江蘇南京 211102)

準(zhǔn)確預(yù)測揮發(fā)性污染物泄漏后在水體和大氣中的時(shí)空分布，能為水氣生態(tài)環(huán)境健康診斷、控制管理提供重要依據(jù)。揮發(fā)性污染物泄漏后在水氣中的時(shí)空分布受到諸多因素影響，如污染物種類，水、氣及其交界面的流態(tài)等。由于水氣界面上揮發(fā)性污染物揮發(fā)傳質(zhì)的復(fù)雜性，迄今為止，人們對(duì)其機(jī)理尚不清楚。本文探討了由風(fēng)速引起的水、氣及其交界面的湍流對(duì)揮發(fā)性污染物揮發(fā)傳質(zhì)的作用。水氣交界面附近的流動(dòng)結(jié)構(gòu)性質(zhì)很獨(dú)特:如分界面通常不平整;在水流速度非常高時(shí)會(huì)形成自然摻氣，水流速度較高時(shí)水面上有波出現(xiàn)，水流速度較低時(shí)，分界面本身變形不大，氣側(cè)的流動(dòng)性質(zhì)與氣流繞過固體邊界的流動(dòng)性質(zhì)相似，而水側(cè)卻并非如此［1-2］。這主要是因?yàn)樗畾饨唤缑嫠畟?cè)的時(shí)均流速梯度很小，因而在水氣交界面附近耗散的湍動(dòng)能只能來源于自由液面以下水側(cè)產(chǎn)生的湍動(dòng)能及水氣交界面以上通過壓強(qiáng)場傳遞的湍動(dòng)能［3］。從這一點(diǎn)來看，水氣交界面附近的流動(dòng)結(jié)構(gòu)極大地依賴于水側(cè)的渦與水氣交界面的碰撞。但當(dāng)氣流速度較高時(shí)，水氣交界面附近的流動(dòng)結(jié)構(gòu)也與氣側(cè)的湍流運(yùn)動(dòng)有關(guān)。

目前還很難通過實(shí)驗(yàn)得到水氣交界面附近流場的準(zhǔn)確數(shù)據(jù)，數(shù)值模擬成為研究水氣交界面湍流結(jié)構(gòu)的有效方法。根據(jù)數(shù)值模擬的精細(xì)程度，湍流數(shù)值模擬可分為直接模擬(DNS)、大渦模擬(LES)和雷諾平均數(shù)值模擬(RANS)［4］。Kawamura等［5-6］將水氣交界面上的速度和剪切應(yīng)力耦合，用DNS方法模擬了水氣交界面的湍流結(jié)構(gòu)，但由于DNS方法對(duì)計(jì)算機(jī)內(nèi)存和計(jì)算速度的要求非常高，目前該方法無法用于工程實(shí)際。Magnaudet等［7］用LES方法研究了水氣交界面上湍流結(jié)構(gòu)及傳質(zhì)過程。Hasegawa等［8-9］提出了DNS和LES相混合的方法:在水氣交界面上的模擬用DNS方法;而水、氣交界面附近的水、氣模擬用 LES方法。Hasegawa等［8-9］雖然考慮了水、氣流動(dòng)的動(dòng)量耦合，但沒有考慮污染物濃度耦合。

本文用氣液兩相流VOF(volume of fluid)方法捕捉水氣交界面［10］，用可實(shí)現(xiàn)k-ε模型計(jì)算水氣交界面湍流結(jié)構(gòu)，在此基礎(chǔ)上利用揮發(fā)性污染物水氣耦合擴(kuò)散模型預(yù)測揮發(fā)性污染物點(diǎn)源泄漏后的時(shí)空分布，分析風(fēng)速引起的水氣湍流對(duì)水氣交界面上揮發(fā)性污染物揮發(fā)傳質(zhì)的作用機(jī)理。

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

1.1 VOF動(dòng)力學(xué)模型

氣液兩相流VOF方法是常見的捕捉水氣交界面的方法，該方法中質(zhì)量守恒、動(dòng)量方程及捕捉水氣交界面體積百分比函數(shù)α方程如下［11］:

式中:ui——網(wǎng)格單元總速度;ρ——網(wǎng)格單元總密度，由液體、氣體的體積加權(quán)平均得出;uLIQi——液體速度;uGASi——?dú)怏w速度;ρGAS——?dú)怏w密度;p——壓強(qiáng);ηt——湍流運(yùn)動(dòng)黏性系數(shù);fj——質(zhì)量力，本文中質(zhì)量力僅有重力。水氣交界面上，0＜α＜1。

1.2 湍流計(jì)算模型

其中

式中:S——應(yīng)變率;k——湍動(dòng)能;ε——湍動(dòng)能耗散率;σm，σn，C2，C1——湍流系數(shù)，σm=1.0，σn=1.2，C2=1.92;η——分子黏性系數(shù)。

ηt由下式確定:

其中

式中:A0，As，U*——系數(shù);Ωi，ji，j——渦通量;ωk—— 角速度;δi，j，k——里奇(Ricci)符號(hào)，j——從角速度為ωk的參考系中觀察到的渦通量。

1.3 揮發(fā)性污染物水氣耦合擴(kuò)散模型

由亨利定律表達(dá)式知，污染物在氣相中的質(zhì)量濃度CGAS與液相中的質(zhì)量濃度CLIQ之比為量綱為1的亨利常數(shù)Haw，即Haw=CGAS/CLIQ。忽略化學(xué)反應(yīng)，揮發(fā)性污染物水氣耦合擴(kuò)散模型表達(dá)式為

式中:C——網(wǎng)格單元總質(zhì)量濃度;Ei，D——污染物在液相、氣相中湍流質(zhì)擴(kuò)散系數(shù)，為ηt與Schmidt數(shù)之比。

2 模擬分析

2.1 模擬條件

計(jì)算域長直水槽14.0 m×0.4 m×0.4 m，空氣部14.0 m×6.4 m×2.6 m。坐標(biāo)原點(diǎn)取在水槽中部靠岸，水槽底部。水流動(dòng)方向?yàn)閤正方向，水槽寬度方向?yàn)閥，垂直向上為z方向。進(jìn)口邊界劃分為水流速度進(jìn)口和空氣速度進(jìn)口，進(jìn)口的k，ε值由經(jīng)驗(yàn)公式給出。假設(shè)出口斷面為充分發(fā)展的湍流，所有物理量一階法向梯度為零，50 mL質(zhì)量濃度為2 g/L的三氯乙烯溶液在x=1.25 m，y=0.20 m處瞬時(shí)泄漏。三氯乙烯的Haw為0.49，在水和空氣中的Schmidt數(shù)分別取754和1.03［13］。在交錯(cuò)網(wǎng)格上，數(shù)值離散用有限體積法，速度壓力耦合用PISO算法，時(shí)間差分采用Crank-Nicholson格式。風(fēng)速太高，水氣交界面將變?yōu)榉沁B續(xù)界面，而本文研究連續(xù)水氣交界面上風(fēng)速引起的湍流對(duì)揮發(fā)性污染物揮發(fā)傳質(zhì)的影響，故模擬風(fēng)速較低，分別取0.4 m/s，0.8 m/s和1.2 m/s。模擬平臺(tái)為開源計(jì)算流體力學(xué)軟件OpenFOAM。

2.2 k和ε與風(fēng)速的關(guān)系

圖1是風(fēng)速為0.8 m/s，1.2 m/s時(shí)，x=3.45 m截面上水氣交界面附近水側(cè)k，ε相對(duì)值隨時(shí)間變化的曲線。如果考察的截面太近，風(fēng)速對(duì)揮發(fā)傳質(zhì)的作用時(shí)間太短，分析無意義。如果考察的截面太遠(yuǎn)，揮發(fā)傳質(zhì)太少，研究意義不大。

圖1 x=3.45 m截面水中 k/kv air=0.4 m/s，ε/εv air=0.4 m/s隨時(shí)間的變化Fig.1 k/kv air=0.4m/s and ε/εv air=0.4m/s as functions of time at cross-section x=3.45 m

圖1中vair=0.4 m/s表示風(fēng)速為0.4 m/s的工況。圖中k/kvair=0.4m/s和ε/εvair=0.4m/s數(shù)值基本大于1，說明風(fēng)速為0.8 m/s，1.2 m/s時(shí)，水側(cè)k，ε大于風(fēng)速0.4 m/s時(shí)水側(cè)k，ε。這是因?yàn)轱L(fēng)速增加，水氣交界面附近時(shí)均速度梯度增加，進(jìn)而導(dǎo)致湍動(dòng)能平均值增加，同時(shí)湍動(dòng)能波動(dòng)幅度提高，氣流誘導(dǎo)的切應(yīng)力及表面變形隨著風(fēng)速增加而促進(jìn)湍能耗散率增加。圖1中(a)，(b)對(duì)應(yīng)的曲線變化趨勢一致，大多數(shù)時(shí)刻，ε/εvair=0.4m/s大于k/kvair=0.4m/s，說明湍能耗散率隨風(fēng)速增加而增加的程度高于湍動(dòng)能隨風(fēng)速增加而增加的程度。

2.3 揮發(fā)性污染物濃度分布

圖2是不同風(fēng)速下，污染物泄漏8 s后，縱截面y=0.225m上污染物質(zhì)量濃度的空間分布。由圖2可見，污染物質(zhì)量濃度分布區(qū)域在縱向上得到了延伸，并且濃度中心由水中經(jīng)過水氣交界面向空氣中轉(zhuǎn)移，進(jìn)而使空氣中的污染物越來越多。水流速度決定污染團(tuán)前鋒的分布狀況，由于3種情況的水流速度相同，故他們的濃度前鋒相差無幾。3種風(fēng)速下污染團(tuán)尾部區(qū)域相差甚遠(yuǎn)。圖2(a)情況下，污染團(tuán)幾乎沒有濃度尾巴;圖2(b)的污染團(tuán)后部拖著長長的濃度尾巴，質(zhì)量濃度值較小;圖2(c)的污染團(tuán)縱向延伸程度最大，濃度尾巴長且數(shù)值高。說明風(fēng)速的增加對(duì)污染團(tuán)尾部的揮發(fā)傳質(zhì)過程影響巨大。由此說明，風(fēng)速引起的k，ε的增加不僅有利于水氣交界面的傳質(zhì)而且有助于污染物在空氣中的擴(kuò)散。這種水氣中污染物濃度相互耦合的特性使得污染團(tuán)尾部的濃度分布占據(jù)的空間隨風(fēng)速增加而增大。傳統(tǒng)非耦合模型認(rèn)為空氣中污染物濃度為零，水中污染物揮發(fā)不受空氣中污染物濃度影響。事實(shí)上，空氣中污染物濃度抑制了水中污染物的揮發(fā)，使得整個(gè)揮發(fā)過程得到延遲。

圖3是x=3.45 m截面水氣交界面附近水側(cè)和氣側(cè)污染物質(zhì)量濃度隨時(shí)間變化的曲線。

圖2 縱截面y=0.225 m污染物質(zhì)量濃度空間分布Fig.2 Pollutant concentration distribution at longitudinal section y=0.225 m

圖3 x=3.45 m截面處不同風(fēng)速下水側(cè)、氣側(cè)污染物質(zhì)量濃度隨時(shí)間的變化Fig.3 Pollutant concentration distribution at cross-section x=3.45 m in water and air at different wind velocities

圖3(a)中，0～4 s期間3條曲線基本沒變化，說明泄漏初期，風(fēng)速對(duì)水中污染物的揮發(fā)傳質(zhì)作用很小。但是從第4秒開始，曲線開始發(fā)生分離，3種風(fēng)速對(duì)應(yīng)的污染物質(zhì)量濃度開始增長，vair=0.8 m/s和vair=1.2 m/s對(duì)應(yīng)的質(zhì)量濃度增長較快。第10秒時(shí)，3種風(fēng)速下污染物質(zhì)量濃度均達(dá)到峰值，vair=1.2 m/s的峰值明顯高于其他2種。在質(zhì)量濃度達(dá)到峰值前，風(fēng)速越高，水中污染物質(zhì)量濃度越高。質(zhì)量濃度達(dá)到峰值后，同一時(shí)刻，vair=0.4 m/s與vair=0.8 m/s對(duì)應(yīng)質(zhì)量濃度幾乎相等?？梢妚air=0.8 m/s時(shí)質(zhì)量濃度下降速率大于vair=0.4 m/s時(shí)質(zhì)量濃度的下降速率。vair=1.2 m/s對(duì)應(yīng)的質(zhì)量濃度稍高些，但其下降速率最大。說明達(dá)到峰值濃度后，風(fēng)速越大，ε越大，水中污染物質(zhì)量濃度下降越快。第6秒時(shí)，圖3(b)中3條曲線開始分離，氣側(cè)污染物質(zhì)量濃度達(dá)到峰值的時(shí)間要早于水側(cè)污染物，說明風(fēng)速對(duì)氣側(cè)污染物擴(kuò)散的影響要先于對(duì)水側(cè)污染物擴(kuò)散的影響。風(fēng)速越大，氣側(cè)污染物質(zhì)量濃度越高，質(zhì)量濃度達(dá)到峰值以后的衰減速率也越大，說明風(fēng)速引起的k，ε的增加有助于水中污染物向空氣中揮發(fā)傳質(zhì)及空氣中污染物的遷移擴(kuò)散。

3 結(jié) 語

將揮發(fā)性污染物在水氣交界面上揮發(fā)傳質(zhì)的耦合特性與可實(shí)現(xiàn)k-ε模型相結(jié)合，揭示湍流特性對(duì)揮發(fā)傳質(zhì)的作用機(jī)理。

風(fēng)速增加，水氣交界面附近時(shí)均速度梯度增加，進(jìn)而導(dǎo)致湍動(dòng)能平均值增加，同時(shí)湍動(dòng)能波動(dòng)幅度提高，氣流誘導(dǎo)的切應(yīng)力及表面變形隨著風(fēng)速增加而促進(jìn)湍能耗散率增加。湍能耗散率隨風(fēng)速增加而增加的程度高于湍動(dòng)能隨風(fēng)速增加而增加的程度。

風(fēng)速引起的湍動(dòng)能和湍能耗散率的提高增強(qiáng)了污染團(tuán)尾部的揮發(fā)傳質(zhì)，風(fēng)速越大，污染團(tuán)尾部的濃度分布占據(jù)的空間越大。風(fēng)速對(duì)氣側(cè)污染物擴(kuò)散的影響要先于對(duì)水側(cè)污染物擴(kuò)散的影響。

［1］LIU S，KERMANI A，SHEN L，et al.Investigation of coupled air-water turbulent boundary layers using direct numerical simulations［J］.Phys Fluids，2009，21:062108.

［2］YAMASHITA S，CHEN C G，TAKAHASHI K，et al.Large scale numerical simulations for multi-phase fluid dynamics with moving interfaces［J］.International Journal of Computation Fluicl Dynamics，2008，22:405-410.

［3］李佳佳，陳春剛，肖鋒.風(fēng)波水氣界面湍流渦旋結(jié)構(gòu)研究［J］.中國科學(xué):物理學(xué) 力學(xué) 天文學(xué)，2011，41:980-994.(LI Jiajia，CHEN Chungang，XIAO Feng.Study on vortex structure near wind-driven air-water interfaces［J］.Sci Sin Phys Mech Astron，2011，41:980-994.(in Chinese))

［4］張兆順，崔桂香，許春曉.湍流理論與模擬［M］.北京:清華大學(xué)出版社，2005:163-165.

［5］KAWAMURA H，OHSAKE K，ABE H，et al.DNS of turbulent heat transfer in channel flow with low to medium-high Prandtl number［J］.International Journal Heat Fluid Flow，1998，19(5):482-491.

［6］LAKEHAL D，F(xiàn)ULGOSI M，YADIGAROGLU G，et al.Direct numerical simulation of turbulent heat transfer across a mobile，sheared gas-liquid interface［J］.ASME Journal Heat Transfer，2003，125(6):1129-1139.

［7］MAGNAUDET J，CALMET I.Turbulent mass transfer through a flat shear-free surface［J］.Fluid Mech Journal，2006，553:155-185.

［8］HASEGAWA Y，KASAGI N.Effects of interfacial velocity boundary condition on turbulent mass transfer at high Schmidt numbers［J］.International Journal Heat Fluid Flow，2007，28(6):1192-1203.

［9］HASEGAWA Y，KASAGI N.Hybrid DNS/LES of high Schmidt number mass transfer across turbulent air-water interface［J］.International Journal of Heat and Mass Transfer，2009，52:1012-1022.

［10］DOLBOW JT，MOSSO S，ROBBINSJ，et al.Coupling volume-of-fluid based interface reconstructions with the extended finite element method ［J］.Computer Methods in Applied Mechanics and Engineering，2008，197:439-447.

［11］CHEN Lipin，JIANG Juncheng.Experiments and numerical simulations on transport of dissolved pollutants around spur dike［J］.Water Science and Engineering，2010，3(3):343-355.

［12］陶文銓.數(shù)值傳熱學(xué)［M］.西安:西安交通大學(xué)出版社，2001:278-280.

［13］施瓦茨巴赫 R P，施格文 P M，英博登 D M.環(huán)境有機(jī)化學(xué)［M］.北京:化學(xué)工業(yè)出版社，2004:425-432.