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

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

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

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

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

1 數(shù)學模型

1.1 VOF動力學模型

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

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

1.2 湍流計算模型

其中

式中:S——應變率;k——湍動能;ε——湍動能耗散率;σ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)符號，j——從角速度為ωk的參考系中觀察到的渦通量。

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

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

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

2 模擬分析

2.1 模擬條件

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

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

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

圖1 x=3.45 m截面水中 k/kv air=0.4 m/s，ε/εv air=0.4 m/s隨時間的變化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表示風速為0.4 m/s的工況。圖中k/kvair=0.4m/s和ε/εvair=0.4m/s數(shù)值基本大于1，說明風速為0.8 m/s，1.2 m/s時，水側(cè)k，ε大于風速0.4 m/s時水側(cè)k，ε。這是因為風速增加，水氣交界面附近時均速度梯度增加，進而導致湍動能平均值增加，同時湍動能波動幅度提高，氣流誘導的切應力及表面變形隨著風速增加而促進湍能耗散率增加。圖1中(a)，(b)對應的曲線變化趨勢一致，大多數(shù)時刻，ε/εvair=0.4m/s大于k/kvair=0.4m/s，說明湍能耗散率隨風速增加而增加的程度高于湍動能隨風速增加而增加的程度。

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

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

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

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

圖3 x=3.45 m截面處不同風速下水側(cè)、氣側(cè)污染物質(zhì)量濃度隨時間的變化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ā)傳質(zhì)作用很小。但是從第4秒開始，曲線開始發(fā)生分離，3種風速對應的污染物質(zhì)量濃度開始增長，vair=0.8 m/s和vair=1.2 m/s對應的質(zhì)量濃度增長較快。第10秒時，3種風速下污染物質(zhì)量濃度均達到峰值，vair=1.2 m/s的峰值明顯高于其他2種。在質(zhì)量濃度達到峰值前，風速越高，水中污染物質(zhì)量濃度越高。質(zhì)量濃度達到峰值后，同一時刻，vair=0.4 m/s與vair=0.8 m/s對應質(zhì)量濃度幾乎相等?？梢妚air=0.8 m/s時質(zhì)量濃度下降速率大于vair=0.4 m/s時質(zhì)量濃度的下降速率。vair=1.2 m/s對應的質(zhì)量濃度稍高些，但其下降速率最大。說明達到峰值濃度后，風速越大，ε越大，水中污染物質(zhì)量濃度下降越快。第6秒時，圖3(b)中3條曲線開始分離，氣側(cè)污染物質(zhì)量濃度達到峰值的時間要早于水側(cè)污染物，說明風速對氣側(cè)污染物擴散的影響要先于對水側(cè)污染物擴散的影響。風速越大，氣側(cè)污染物質(zhì)量濃度越高，質(zhì)量濃度達到峰值以后的衰減速率也越大，說明風速引起的k，ε的增加有助于水中污染物向空氣中揮發(fā)傳質(zhì)及空氣中污染物的遷移擴散。

3 結(jié) 語

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

風速增加，水氣交界面附近時均速度梯度增加，進而導致湍動能平均值增加，同時湍動能波動幅度提高，氣流誘導的切應力及表面變形隨著風速增加而促進湍能耗散率增加。湍能耗散率隨風速增加而增加的程度高于湍動能隨風速增加而增加的程度。

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

［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］李佳佳，陳春剛，肖鋒.風波水氣界面湍流渦旋結(jié)構(gòu)研究［J］.中國科學:物理學 力學 天文學，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］.北京:清華大學出版社，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ù)值傳熱學［M］.西安:西安交通大學出版社，2001:278-280.

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