陈丽萍，程璟涛，蒋军成，邓广发

(1.南京工业大学城市建设与安全工程学院，江苏南京 210009;2.江苏方天电力技术有限公司，江苏南京 211102)

准确预测挥发性污染物泄漏后在水体和大气中的时空分布，能为水气生态环境健康诊断、控制管理提供重要依据。挥发性污染物泄漏后在水气中的时空分布受到诸多因素影响，如污染物种类，水、气及其交界面的流态等。由于水气界面上挥发性污染物挥发传质的复杂性，迄今为止，人们对其机理尚不清楚。本文探讨了由风速引起的水、气及其交界面的湍流对挥发性污染物挥发传质的作用。水气交界面附近的流动结构性质很独特:如分界面通常不平整;在水流速度非常高时会形成自然掺气，水流速度较高时水面上有波出现，水流速度较低时，分界面本身变形不大，气侧的流动性质与气流绕过固体边界的流动性质相似，而水侧却并非如此［1-2］。这主要是因为水气交界面水侧的时均流速梯度很小，因而在水气交界面附近耗散的湍动能只能来源于自由液面以下水侧产生的湍动能及水气交界面以上通过压强场传递的湍动能［3］。从这一点来看，水气交界面附近的流动结构极大地依赖于水侧的涡与水气交界面的碰撞。但当气流速度较高时，水气交界面附近的流动结构也与气侧的湍流运动有关。

目前还很难通过实验得到水气交界面附近流场的准确数据，数值模拟成为研究水气交界面湍流结构的有效方法。根据数值模拟的精细程度，湍流数值模拟可分为直接模拟(DNS)、大涡模拟(LES)和雷诺平均数值模拟(RANS)［4］。Kawamura等［5-6］将水气交界面上的速度和剪切应力耦合，用DNS方法模拟了水气交界面的湍流结构，但由于DNS方法对计算机内存和计算速度的要求非常高，目前该方法无法用于工程实际。Magnaudet等［7］用LES方法研究了水气交界面上湍流结构及传质过程。Hasegawa等［8-9］提出了DNS和LES相混合的方法:在水气交界面上的模拟用DNS方法;而水、气交界面附近的水、气模拟用 LES方法。Hasegawa等［8-9］虽然考虑了水、气流动的动量耦合，但没有考虑污染物浓度耦合。

本文用气液两相流VOF(volume of fluid)方法捕捉水气交界面［10］，用可实现k-ε模型计算水气交界面湍流结构，在此基础上利用挥发性污染物水气耦合扩散模型预测挥发性污染物点源泄漏后的时空分布，分析风速引起的水气湍流对水气交界面上挥发性污染物挥发传质的作用机理。

1 数学模型

1.1 VOF动力学模型

气液两相流VOF方法是常见的捕捉水气交界面的方法，该方法中质量守恒、动量方程及捕捉水气交界面体积百分比函数α方程如下［11］:

式中:ui——网格单元总速度;ρ——网格单元总密度，由液体、气体的体积加权平均得出;uLIQi——液体速度;uGASi——气体速度;ρGAS——气体密度;p——压强;ηt——湍流运动黏性系数;fj——质量力，本文中质量力仅有重力。水气交界面上，0＜α＜1。

1.2 湍流计算模型

其中

式中:S——应变率;k——湍动能;ε——湍动能耗散率;σm，σn，C2，C1——湍流系数，σm=1.0，σn=1.2，C2=1.92;η——分子黏性系数。

ηt由下式确定:

其中

式中:A0，As，U*——系数;Ωi，ji，j——涡通量;ωk—— 角速度;δi，j，k——里奇(Ricci)符号，j——从角速度为ωk的参考系中观察到的涡通量。

1.3 挥发性污染物水气耦合扩散模型

由亨利定律表达式知，污染物在气相中的质量浓度CGAS与液相中的质量浓度CLIQ之比为量纲为1的亨利常数Haw，即Haw=CGAS/CLIQ。忽略化学反应，挥发性污染物水气耦合扩散模型表达式为

式中:C——网格单元总质量浓度;Ei，D——污染物在液相、气相中湍流质扩散系数，为ηt与Schmidt数之比。

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，ε值由经验公式给出。假设出口断面为充分发展的湍流，所有物理量一阶法向梯度为零，50 mL质量浓度为2 g/L的三氯乙烯溶液在x=1.25 m，y=0.20 m处瞬时泄漏。三氯乙烯的Haw为0.49，在水和空气中的Schmidt数分别取754和1.03［13］。在交错网格上，数值离散用有限体积法，速度压力耦合用PISO算法，时间差分采用Crank-Nicholson格式。风速太高，水气交界面将变为非连续界面，而本文研究连续水气交界面上风速引起的湍流对挥发性污染物挥发传质的影响，故模拟风速较低，分别取0.4 m/s，0.8 m/s和1.2 m/s。模拟平台为开源计算流体力学软件OpenFOAM。

2.2 k和ε与风速的关系

图1是风速为0.8 m/s，1.2 m/s时，x=3.45 m截面上水气交界面附近水侧k，ε相对值随时间变化的曲线。如果考察的截面太近，风速对挥发传质的作用时间太短，分析无意义。如果考察的截面太远，挥发传质太少，研究意义不大。

图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数值基本大于1，说明风速为0.8 m/s，1.2 m/s时，水侧k，ε大于风速0.4 m/s时水侧k，ε。这是因为风速增加，水气交界面附近时均速度梯度增加，进而导致湍动能平均值增加，同时湍动能波动幅度提高，气流诱导的切应力及表面变形随着风速增加而促进湍能耗散率增加。图1中(a)，(b)对应的曲线变化趋势一致，大多数时刻，ε/εvair=0.4m/s大于k/kvair=0.4m/s，说明湍能耗散率随风速增加而增加的程度高于湍动能随风速增加而增加的程度。

2.3 挥发性污染物浓度分布

图2是不同风速下，污染物泄漏8 s后，纵截面y=0.225m上污染物质量浓度的空间分布。由图2可见，污染物质量浓度分布区域在纵向上得到了延伸，并且浓度中心由水中经过水气交界面向空气中转移，进而使空气中的污染物越来越多。水流速度决定污染团前锋的分布状况，由于3种情况的水流速度相同，故他们的浓度前锋相差无几。3种风速下污染团尾部区域相差甚远。图2(a)情况下，污染团几乎没有浓度尾巴;图2(b)的污染团后部拖着长长的浓度尾巴，质量浓度值较小;图2(c)的污染团纵向延伸程度最大，浓度尾巴长且数值高。说明风速的增加对污染团尾部的挥发传质过程影响巨大。由此说明，风速引起的k，ε的增加不仅有利于水气交界面的传质而且有助于污染物在空气中的扩散。这种水气中污染物浓度相互耦合的特性使得污染团尾部的浓度分布占据的空间随风速增加而增大。传统非耦合模型认为空气中污染物浓度为零，水中污染物挥发不受空气中污染物浓度影响。事实上，空气中污染物浓度抑制了水中污染物的挥发，使得整个挥发过程得到延迟。

图3是x=3.45 m截面水气交界面附近水侧和气侧污染物质量浓度随时间变化的曲线。

图2 纵截面y=0.225 m污染物质量浓度空间分布Fig.2 Pollutant concentration distribution at longitudinal section y=0.225 m

图3 x=3.45 m截面处不同风速下水侧、气侧污染物质量浓度随时间的变化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条曲线基本没变化，说明泄漏初期，风速对水中污染物的挥发传质作用很小。但是从第4秒开始，曲线开始发生分离，3种风速对应的污染物质量浓度开始增长，vair=0.8 m/s和vair=1.2 m/s对应的质量浓度增长较快。第10秒时，3种风速下污染物质量浓度均达到峰值，vair=1.2 m/s的峰值明显高于其他2种。在质量浓度达到峰值前，风速越高，水中污染物质量浓度越高。质量浓度达到峰值后，同一时刻，vair=0.4 m/s与vair=0.8 m/s对应质量浓度几乎相等。可见vair=0.8 m/s时质量浓度下降速率大于vair=0.4 m/s时质量浓度的下降速率。vair=1.2 m/s对应的质量浓度稍高些，但其下降速率最大。说明达到峰值浓度后，风速越大，ε越大，水中污染物质量浓度下降越快。第6秒时，图3(b)中3条曲线开始分离，气侧污染物质量浓度达到峰值的时间要早于水侧污染物，说明风速对气侧污染物扩散的影响要先于对水侧污染物扩散的影响。风速越大，气侧污染物质量浓度越高，质量浓度达到峰值以后的衰减速率也越大，说明风速引起的k，ε的增加有助于水中污染物向空气中挥发传质及空气中污染物的迁移扩散。

3 结 语

将挥发性污染物在水气交界面上挥发传质的耦合特性与可实现k-ε模型相结合，揭示湍流特性对挥发传质的作用机理。

风速增加，水气交界面附近时均速度梯度增加，进而导致湍动能平均值增加，同时湍动能波动幅度提高，气流诱导的切应力及表面变形随着风速增加而促进湍能耗散率增加。湍能耗散率随风速增加而增加的程度高于湍动能随风速增加而增加的程度。

风速引起的湍动能和湍能耗散率的提高增强了污染团尾部的挥发传质，风速越大，污染团尾部的浓度分布占据的空间越大。风速对气侧污染物扩散的影响要先于对水侧污染物扩散的影响。

［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］李佳佳，陈春刚，肖锋.风波水气界面湍流涡旋结构研究［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，FULGOSI 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］陶文铨.数值传热学［M］.西安:西安交通大学出版社，2001:278-280.

［13］施瓦茨巴赫 R P，施格文 P M，英博登 D M.环境有机化学［M］.北京:化学工业出版社，2004:425-432.