90°等宽明渠交汇水流污染物混合规律数值模拟

毛颂平，朱 海，王玲玲，华祖林，陆圣杰，薛晓鹏

(1.河海大学水利水电学院，江苏 南京 210098； 2.河海大学环境学院，江苏 南京 210098；3.宁波市水利水电规划设计研究院有限公司，浙江 宁波 315192)

交汇水流是天然河流及人工河道中普遍存在的水力现象。干支流交汇处水流紊动剧烈，对污染物的混合特性有着重要影响。随着水环境问题的日益突出，国内外学者对明渠交汇水流污染物混合规律进行了大量研究。在物理试验研究方面，陈凯霖等[1]基于非等宽交汇口模型研究了污染物质量浓度场的分布规律，探究了交汇角和污染物质量浓度对其的影响；Ribeiro等[2]采用野外实测和水槽模型试验相结合的方式，研究发现干支流流量比对交汇口物质运输有着重要影响，在流量比不变的情况下，拓宽支流宽度有利于污染物混合；Hua等[3]通过建立辫状多支流河流物理模型研究污染物的混合特性，研究表明污染物的释放位置是污染物迁移过程的主要影响因素；唐洪武等[4]通过实测发现悬沙质量浓度与水温更能反映长江与鄱阳湖交汇前后的水质掺混特性。在数值模拟研究方面，刘雪兰[5]应用SMA二维模型模拟了重庆两江段流场和水质，发现在汇流口一侧的河岸会形成高质量浓度的污染带；李娟等[6]建立二维水动力水质模型模拟汾河入黄口污染带分布特性，研究表明横向上扩散起主导作用，纵向上对流起主导作用；米潭等[7]建立平面二维自由表面流MIKE 21 FM模型，研究了Y形交汇区污染物扩散规律，得出了汇流比与污染物横向、纵向质量浓度梯度的关系；Lee等[8]建立二位对流弥散数学模型，研究了韩国汉江某分汊河道下游的水面形态与污染物输运过程；茅泽育等[9]建立了明渠交汇口水流流动及污染物输移模型，较好地模拟了污染混合区的大小；Biron等[10]分别采用试验模型和天然河流模拟了交汇口污染物的混合，发现干支渠底部存在高差会加速污染物的混合；魏娟等[11]建立水气两相流数学模型，研究了不同汇流比工况下交汇口污染物输移扩散特性，得出汇流比越大，污染带越狭长的结论；Ramon等[12]通过三维数学模型分析了埃布罗河和塞格雷河交汇口污染物的混合速率；Tang等[13]采用三维数值模型研究发现污染物的混合受水流形态的强烈影响，而汇流比和河床形态又深刻影响着水流形态；顾莉等[14]通过建立U形交汇河道水气两相流模型研究发现污染物横向离散系数沿程呈单峰结构分布，而纵向离散系数呈双峰结构分布。

综上可知，汇流比是影响交汇水流水动力特性的重要因素，对交汇河口污染物输运及混合特性具有重要影响，目前针对90°等宽明渠交汇水流污染带三维空间分布特征的定量研究仍有不足，因此本文基于三维水动力-污染物耦合输运数学模型，着重研究了不同汇流比情况下明渠交汇口下游污染物三维空间分布特征，分析了污染物质量浓度的混合特性，获得了污染带最大宽度与汇流比的定量关系。

1 数学模型

1.1 控制方程

本文采用基于雷诺时均方程(RANS)的三维水气两相流模型，模拟研究区域污染物质量浓度的分布。计算区域内水流运动满足质量守恒方程及动量守恒方程，污染物质量浓度满足质量浓度守恒方程。本文采用雷诺应力模型封闭控制质量守恒方程和动量守恒方程，根据时均化法则，直接建立并模化雷诺应力的输运方程，没有采用涡黏性假设，可以有效地计算各向异性的湍流流场，对复杂流动的模拟有着更强的优势。

雷诺应力输运方程如下：

DT,ij+DL,ij+Pij+Πij+εij

(1)

该方程需联合耗散率的输运方程进行求解：

(2)

本文采用质量浓度通量的二阶矩模式封闭控制质量浓度守恒方程：

(3)

1.2 数值方法及边界条件

控制方程的离散方法采用有限体积法，控制方程组扩散项的离散均采用中心差分格式，水动力模块控制方程对流项的离散采用QUICK格式，污染物输运模块控制方程对流项的离散采用二阶迎风格式。采用PISO算法对压力和速度场进行耦合计算。

交汇水流的相互作用导致交汇口水面剧烈变化，而自由表面的模拟对于模型的可靠性与准确性十分重要，本文采用VOF法模拟自由表面。计算区域底部及侧壁采用无滑移边界条件。主、支渠入口边界采用速度进口，流向垂直于进口断面。出口边界采用压力出口，下游出口水位保持与试验值一致，出口压力满足静水压强分布，液面相对压强为0。进出口边界处紊动能和耗散率的值均根据经验公式[11]确定：k=0.003 75u2，ε=k1.5/(0.4H)，其中u为断面平均速度，H为水深。

2 模型验证

2.1 水动力模型验证

本文将水动力模型数值模拟结果与Weber等[17]的90°等宽明渠交汇流物理试验数据作对比，验证工况下游出口尾水高度H0= 0.306 m，弗劳德数Fr=0.37，雷诺数Re=186 000，汇流比q=Qu/Qd= 0.25(Qu为主渠进口流量，Qd为下游出口流量，Qd= 0.17 m3/s)。数学模型计算区域为矩形断面明渠交汇口，主、支渠交汇角为90°，沿程等宽，底坡水平，采用顺接方式连接，主、支渠适当加长是为了消除边界条件的影响，使水流紊动充分发展，计算区域如图1所示。采用结构化六面体网格剖分计算区域，对交汇口、近壁区域、水气交界面和渠底网格进行加密。经过网格无关性分析，验证工况网格数为24.6万。

图1 计算区域示意图

数学模型弗劳德数、雷诺数及流量比均与物理试验保持一致，坐标系以渠宽W进行无量纲化处理，速度以下游出口平均流速U0进行无量纲化处理。水面线模拟值与试验值对比结果表明，不同纵向断面水面线高度与物理试验实测值吻合较好，交汇口附近不同测点的速度数值模拟值与物理试验实测值吻合良好，能够反映流速变化的整体趋势。经过验证，该模型能够较好地反映交汇口水流的水动力特性。典型纵断面水面线对比如图2所示，典型断面流速对比如图3所示。

图2 y/W=0.5纵断面水面线对比

图3 x/W=-1横断面各测线上的流速对比

2.2 污染物输运模型验证

为了验证数值模型对于污染物输运的可靠性，本文选取陈凯霖等[1,18]的物理试验数据进行对比。研究对象主、支渠交汇角为90°，主渠宽度为W，支渠宽度为1/3W，主、支渠上游长10W，主渠下游长10W。计算区域采用六面体网格进行剖分，总网格数为19.2万。计算工况Fr=0.10，Re=6 357.79，q=0.81，主流入口污染物质量浓度为0 μg/L，支流入口污染物质量浓度为2 000 μg/L，均与物理试验保持一致。图4为验证工况下x/W=0.8和x/W=1.33横断面上y/W=0.5测线上流速实测值与计算值的对比，可见流速模拟值与测量值分布趋势基本一致，由于本文采用的是雷诺时均模型，模拟获得的流速分布曲线相对试验值更为光滑平顺，总体误差在10%以内。图5为验证工况下y/W=0.167和y/W=0.5纵断面上z/W=0.427深度处污染物质量浓度纵向分布计算值与实测值的对比，两者分布趋势吻合基本良好。

图4 不同横断面y/W=0.5处流速对比

图5 不同纵断面z/W=0.427处污染物质量浓度对比

3 不同汇流比条件下明渠交汇水流污染物三维分布规律

3.1 工况设置

为研究明渠交汇水流污染物质量浓度的三维分布规律及汇流比对其的影响，采用验证后的模型对q=0.083、0.25、0.417、0.583、0.750、0.917等6种工况进行计算，所选取的工况与Weber等[17]的90°等宽明渠交汇流物理试验选取工况一致。所有工况主、支渠等宽，交汇角均为90°，下游出口断面Fr=0.37，Re=186 000，Qd=0.17 m3/s。主流入口污染物质量浓度C=0 μg/L，支流入口污染物质量浓度C0=1 μg/L。

3.2 污染物质量浓度三维空间分布特征

图6和图7分别展示了汇流比为0.25和0.75时支流污染物进入主渠后下游各横断面的质量浓度分布。支流高质量浓度污染物进入主渠后，在主渠靠近支渠一侧形成污染物高质量浓度区域；污染物随着支流的初始动能冲向主渠右岸，沿主渠横向断面形成明显的质量浓度梯度；交汇口下游x/W>-3范围内污染物等质量浓度线总体保持垂直，在垂向上质量浓度变化不明显；随着干支流的进一步混合，远离交汇口的下游区域(x/W<-3)等质量浓度线逐渐发生倾斜和扭曲，污染物质量浓度呈现明显的三维特征。

图6 q=0.25时交汇口下游典型横断面流速矢量及污染物质量浓度分布云图

图7 q=0.75时交汇口下游典型横断面流速矢量及污染物质量浓度分布云图

模拟结果表明，不同汇流比条件下交汇口下游污染物空间分布特征具有显著差异。不同汇流比条件下污染物横向输运范围不同，汇流比较小时，污染物横向输运的范围更大，汇流比较大时，污染物横向输运的范围较小，这主要是不同汇流比工况下支流进入主渠的动能不同导致的；不同水深处污染物质量浓度横向混合范围也和汇流比相关，汇流比较小时，表层污染物横向混合范围大于底层污染物混合范围(图6)，污染物质量浓度等值线呈D形分布；汇流比较大时，表层污染物横向混合范围小于底层污染物混合范围(图7)，污染物质量浓度等值线呈L形分布。污染物质量浓度分布特征与断面环流紧密相关，在小汇流比工况下，由于在渠道中部存在顺时针的断面环流(图6(c)(d))，靠近渠道左岸的污染物向表层水体输运，靠近渠道右岸的污染物向底层水体输运，等质量浓度线在水体表层和底层向左岸移动，在水体中部向右岸移动；大汇流比工况下，由于在靠近渠道左岸处存在逆时针的断面环流(图7(c)(d))，靠近渠道左岸的污染物向底层水体输运，等质量浓度线逐渐发生倾斜，逐渐形成L形的高质量浓度区域。

3.3 污染物质量浓度横向分布范围

污染带的最大宽度表征着污染物横向分布的最大范围，对于研究污染物对水体的影响范围有着重要意义。在一般情况下，当边缘点的质量浓度为同一断面上最大质量浓度的5%时，该点被认为是污染带的边界点[19]，因此污染带最大宽度b定义为污染带边缘质量浓度是断面质量浓度5%时所占有的最大宽度(如图8(d)所示)。图8为不同汇流比工况表层水体(z/W=0.3)污染物质量浓度分布云图，图中红色虚线为污染带边界。当支流流量占比较大时(对应q≤0.417)，表层水体污染带总能在交汇口下游x/W=-2范围内影响到右岸(图8(a)(b)(c))；对于支流流量占比较小时(q>0.417)，污染带最大宽度b随着汇流比的增大呈逐渐减小趋势(图8(d)(e)(f))。

图8 不同汇流比条件下交汇口表层污染物质量浓度分布云图

数值计算结果表明，不同水深处污染带宽度与汇流比具有不同的相关关系。图9为各工况下不同水深处污染带相对最大宽度b/W与汇流比的关系曲线。由图9(a)可知，交汇口表层水体(z/W=0.3)的污染物在小汇流比工况下(如q<0.417)其横向污染范围已达到右岸，污染带最大宽度不随汇流比的变化而变化；在大汇流比工况(如q>0.417)，交汇口表层水体污染带最大宽度随着汇流比的增大而减小，呈现良好的一次函数关系。中下层水体污染带最大宽度分布与表层水体不同，统计不同水深处各工况污染带最大宽度，并作无量纲化处理，拟合的曲线表明b/W与汇流比在中下层水体中呈现相似的规律且有着良好的二次曲线关系，见图9(b)。

图9 各工况下不同水深污染带b/W与q的关系曲线

3.4 污染物质量浓度混合特性

通常用断面污染物质量浓度的不均匀系数δ及标准差σ的沿程分布评估污染物的混合速率[12]。对于一个给定的断面，δ表征与下游完全混合后污染物质量浓度的偏差，计算公式如下：

δ=(Ci-Cp)/Cp

(4)

式中：Ci为某断面污染物平均质量浓度计算值；Cp为污染物平均预测质量浓度。δ越接近于0，污染物混合越充分。污染物的平均预测质量浓度Cp，即充分混合之后河道内污染物的质量浓度可按下式计算：

Cp=(CmQm+CtQt)/(Qm+Qt)

(5)

式中：Cm和Ct分别为主流和支流的污染物质量浓度；Qm和Qt分别为主流和支流的流量。

图10为不同汇流比工况下交汇口下游各横断面δ的沿程分布。由图10可见，各工况污染物不均匀系数的沿程变化呈现相似的变化规律：先增大，在x/W=-2断面达到最大值后逐渐减小，减小的速率先快后慢，最后趋于稳定。汇流比越大，δ的峰值越大，即断面平均质量浓度与平均预测质量浓度的差距越大；汇流比越大，δ的下降速度也越大，即污染物混合的速率越快。

图10 交汇口下游各横断面δ变化曲线

δ反映的是给定断面污染物质量浓度与完全混合后水体中污染物质量浓度的偏差，而给定断面污染物质量浓度分布均匀程度通常用标准差σ来表征[12]。图11为不同工况下交汇口下游主渠各横断面标准差σ的变化曲线。由图11可知，干支流流量越接近(汇流比越接近0.5)，标准差越大，污染物混合越不充分。在大汇流比情况下(如q=0.583、0.750、0.917)，标准差呈现先迅速增大后缓慢减小的规律；在小汇流比工况下(如q=0.083、0.250、0.417)，标准差呈现先迅速增大后减小，再小幅增大，最后缓慢减小的规律。这是由于小汇流比工况下，主渠近支渠一侧形成了明显的分离区[20]。分离区的存在使得过流面积收缩，主流在经过分离区时流速先增大后减小。因此，在x/W=-2断面下游由于流速缓慢减小，污染物向下游输运的质量浓度通量有所降低，从而减弱了污染物的混合程度，故标准差出现小幅增大。

图11 交汇口下游各横断面σ变化曲线

4 结 论

a.交汇水流污染物分布三维特征明显，汇流比对交汇水流污染物的分布有着重要影响，主要体现在污染物横向分布上。汇流比较小时，污染物横向输运的范围较大，汇流比较大时，污染物横向输运的范围较小；随着干支流的混合，交汇口下游断面等质量浓度线逐渐发生倾斜与扭曲。

b.污染物的分布与二次流密切相关。小汇流比工况下，渠道中部存在顺时针的断面环流，交汇口下游各断面污染物质量浓度总体呈中部质量浓度高、表层及底层质量浓度低的D形分布；大汇流比工况下，渠道中部存在逆时针的断面环流，交汇口下游各断面污染物质量浓度总体呈底层质量浓度高、中上层质量浓度低的L形分布。

c.小汇流比工况下，表层水体污染带最大宽度为主渠宽度，不随汇流比的变化而改变，大汇流比工况下，污染带最大宽度与汇流比呈现良好的一次函数关系；中下层水体污染带最大宽度随着汇流比的增大而减小，与汇流比呈二次曲线关系。

d.污染物的混合速率与混合程度均与汇流比密切相关。汇流比越大，污染物混合的速率越快，汇流比越接近0.5，污染物越不易充分混合；小汇流比情景下交汇口下游出现的大尺度分离区会在局部区域内增加污染物质量浓度分布的不均匀性。