耦合VOF法的河道三维水流特性数值模拟研究

周正印，杨 楠

（天津市市政工程设计研究院，天津市 300392）

1 研究背景

随着我国经济的深入发展，国家对生态环境的保护越来越重视，污染防治是2018年我国提出的三大攻坚战之一，被放在了前所未有的重要位置上。但是长期以来，我国污染治理进程相对滞后，河湖空间管理和控制能力相对薄弱，水流不畅、水岸杂乱、水景不佳、水体黑臭等水环境问题日趋恶化，内涝等水安全问题日益凸显，成为经济社会和生态环境协调发展的制约因素。城市中小河道水体黑臭是点源污染、面源污染、内源污染等综合因素导致的结果[1]，其中丧失生态功能的水体，往往流动性降低或完全消失，导致水体复氧能力衰退，形成缺氧、厌氧的水环境，引发水质恶化，加剧水体黑臭。目前针对黑臭水体的治理方法主要是包括截污、清淤、生态、曝气、活水等一系列方法在内的综合处理[2]，而其中河道水动力条件是实施上述治理措施的依据，因此对河道水流进行数值模拟，改善河道水动力条件，可为黑臭水体治理提供前期支持。

目前，对河道水流特性进行数值模拟的研究起步较早，研究较多。在国外，Renato[3]采用二维高精度浅水方程对意大利2014年Secchia河的洪水事件进行了数值模拟，为达到较快的计算速度，采用了并行计算方法；Siavash[4]以伊朗Karun河的实测数据为基础，采用二维数值模型，将固定的曼宁糙率改为随河道条件相应动态变化的糙率，提高了河道水流数值模拟的精度；Ao[5]考虑水质，采用三维数值模型对河道水流进行了数值模拟，为河道水质管理提供了技术支持。在国内，张玮[6]针对目前大多河流存在着平面顺直、断面单一的问题，以奥地利某顺直河道段低水生态修复工程为例，利用二维水流数学模型，研究修复工程的平面布局；彭毕帅[7]采用二维浅水方程对长江下游黑沙洲水道在不同流量下滩脊线加高、岸线崩退、左槽封堵3中工况下的水流变化进行了数值模拟；王甲荣[8]建立基于有限体积法的平面二维水动力模型对海河流域鲁北平原典型河流弯道进行数值模拟，分析不同工况下河床、河道边坡的冲淤变化，选出适宜河流弯道生态治理的河床形态。

综上所述，目前国内外对河流水流特性的数值模型多集中于一、二维，忽略了对河流自由水面问题的研究，并对水流竖直方向上的变化进行了平均化处理，缺乏对河流自由水面波动的模拟，简化了实际河流的复杂弯曲流态与水流细节，降低了模拟的精度。针对上述问题，本文建立了三维高精度河道地形模型，考虑河道水流的三维非稳态运动，建立耦合流体体积函数（volume of fluid，VOF）法与紊流模型的三维数值模型，研究了河道干流段水流三维运动规律，得到各断面水流的水深、流速等特征值，为进一步改善河流水动力条件及水生态措施的实施提供理论依据和决策支持。

2 研究方法

2.1 耦合VOF法的三维k-ε模型

流体力学中浅水流动模型（即圣维南方程组）是对实际流动的一种简化和概化的数学模型，无法精确、稳定地模拟河道水流问题[9]。因此，本文提出一种耦合VOF法和k-ε双方程紊流模型的数值模拟方法，使之能够适用于求解更广泛的具有复杂边界形态和流动特征的自由表面流动。

VOF法允许较陡的自由表面和非单一表面，对河流水质迁移转化自由水面数值模拟采用VOF法可精确计算各计算断面水位、自由水面形状、污染物分布及流速等变量，并以可视化形式表达。VOF法通过定义一个流体体积函数F，F是空间和时间的函数，即 F（x，y，z，t）。在离散网格内，F 取值是网格内各相流体的体积与能够被流体通过的空间体积的比值。在每一个单元内，变量或代表一相，或代表多相，用aq表示单元内第q相流体的体积分数。在每一个控制体中，所有的相体积分数之和为1，即：

对于第q相流体，体积分数连续方程为：

式中：u、v、w 分别为主场沿 x、y、z方向的分速度；Saq为该相的质量源，无源情况下为零。

k-ε紊流模型的连续性方程、动量方程、紊动能k方程和ε方程如下：

连续性方程：

式中：t为时间，s；ui为速度分量，m/s；xi为坐标分量；ρ为密度，kg/m3；μ 为分子粘性系数，N·m/s；P为修正压力，Pa；ui为紊流粘性系数；CD为源项；k为紊动动能，m2/s2；ε 为紊动耗散率，m2/s2；∂k、∂ε 分别为 k、ε 的紊流普朗特数，无因次；C1ε、C2ε为经验常数，无因次[10]。

2.2 模型的边界条件

（1）进口边界条件：如果进口条件已知，给定速度如下式，湍流参数和根据经验公式计算得到；如果不知道准确的来流条件，只知道来流的流量条件（Q为进口流量），A为进口面积，可以速度方向给定在结构化网格方向上（如河道），并设m→为网格的方向矢量，则进口速度给定如下：

（2）出口边界条件：通常在计算域的出口，各速度分量(u，v，w)以及k和ε均取为第二类齐次边界条件，即（n→代表出口断面的法向）：

（3）固壁边界条件：在壁面上采用无滑移条件，即：u=v=w=0。为计算近壁区的紊流，采用壁面函数法。

（4）压力边界条件：在计算域的边界上，压力应满足Neumann条件，即第二类边界条件[11]。为保证计算的稳定性，在规定的某一内点上，压力为给定值，且定义为参考压力Pref。

2.3 模型求解

采用有限体积法离散控制方程，差分格式采用QUICK差分格式[12]，并采用基于SIMPLE算法改进的压力的隐式算子分割算法（Pressure Implicit Split Operator，PISO 算法）对方程进行求解[13]。

3 模型验证

针对三维模型的特点，采用S.Soare Frazao与Y.Zech在急转弯河道模型中进行的试验结果作为验证的依据[14]。

试验在一个有90°急弯的模型河道中进行。上游设置一个水箱用于储水，平面大小为2.44 m×2.44 m×2.39 m，河道为矩形断面，宽0.495 m，玻璃壁面；上游河道长3.92 m，急弯下游河道长2.92 m，河道床面高出水库底面0.33 m。末端自由出流，水箱中的初始水位高出河道床面0.25 m，实验布置见图1。模型采用贴体网格技术对模型进行划分，网格单元的大小为x×y×z=0.01 m×0.01 m×0.01 m，共划分出88 130个网格，计算网格见图2。

图1 溃坝实验模型示意图(单位：m)

图2 网格模型

实验前河道为干河床，试验时将水箱闸门快速提起，模拟河水在河道中的流动，通过玻璃边壁处的摄影捕捉某些断面的水位变化过程，图3为沿程水深模拟值与实验值对比。将模拟得到的t=3 s与t=7 s时下有河道中水深与实验值对比，模拟结果与实验结果吻合，最大误差为9.7%，最小误差为3.1%，平均误差为6.5%。

4 工程应用

4.1 工程简介

阜阳市位于安徽省西北部，淮河以北，是淮海经济区重要组成部分。近年来，阜阳市的经济建设取得了长足的发展，但是污染治理进程相对滞后，由于合流制管道排污、初期雨水污染、水体流动性差、污染物淤积等因素，导致部分河道水体黑臭。为改善阜阳市城区人居环境，提升城市品位，阜阳市决定启动阜阳市城区水系综合整治（含黑臭水体治理）项目。本文以其中某典型黑臭水体河道为研究对象，应用建立的耦合VOF法和双方程的紊流数值模型对河道水流进行数值模拟，为进一步的黑臭水体综合治理措施提供科学依据和决策支持。

图3 沿程水深模拟值与实验值对比(单位：m)

4.2 计算模型及网格划分

该河道长5.2 km，河水自西向东流动，上游为盲端，下游接入另一条河道，中间与3条城区河道交叉，加上沿河雨水及合流制排口，为该河道的主要来水。依照河道平面CAD布置图及断面图，将包含河底形状及高程等实际地形数据的格式文件导入到CFD软件中，实现了真实地形资料在CFD软件计算网格模型中的精确表达，弥补了以往CFD建模中通过坐标绘制网格模型而使网格精确性不足的局限。采用三角形网格、无结构贴体网格划分技术和局部加密网格生成技术建立干流河道的计算网格模型，河道网格单元总数为712718个。河道流域及河道干流的三维网格模型的全局图见图4。由于河道长宽高比例大，因此在全局显示的情况下将使得河道在Y、Z方向网格几乎无法清楚显示，河道几乎成为线状，但河道各弯道及转弯仍可见。

根据该河道设计要求，确定该河道设计流量为10.40 m3，河道一设计流量为11.20 m3，河道二设计流量为22.10 m3，河道三设计流量为30.10 m3，水动力计算采用进口流速边界，出口采用流速边界，并保证总进口流量等于总出口流量，河道基本达到稳定状态。

图4 干流河道网格模型

为了能够更好的分析河道水力条件，在河道干流设置了8个典型断面监测水流变量，分别在河道起始点、河道交叉处、河道中间、河道转弯、河道末端等处，断面详细布置情况见图5。

图5 河道干流8个典型断面分布

4.3 VOF动态及水流流速分析

图6为干流河道已基本达到稳定状态时，6个典型断面的VOF分布情况。根据图3模拟结果计算水深可得，E1、E2、E3、E4、E5、E6、E7 和 E8 断面的平均水深值分别为2.79 m、2.98 m、3.00 m、3.02 m、3.30 m、3.36 m、3.42 m和3.43 m，断面处的水深分布见图7。

由图7可以看出，随着断面从上游向下游推进，河道水深逐渐变深，并且断面水平面的波动比上游剧烈。断面1由于处于河道起点，且为盲端，没有其他补充水源，水平面较为平稳，EI断面最小水深值为2.790 m，最大水深值为2.791 m；E2断面位于本研究河道与河道1交叉处，河道断面具有一定的波动性，水深比E1断面大，最小水深值为2.978 m，最大水深值为2.983 m；E3断面位于本研究河道与河道2交叉处，水深比E2断面大，最小水深值为2.997 m，最大水深值为3.003 m；E4断面位于河道顺直段，水流较为平缓，水位变化不大，最小水深值为3.019 m，最大水深值为3.020 m；E5断面位于本研究河道与和河道3交叉口处，且该处河道转弯，水流变化最为剧烈，水深增加较大，最小水深值为3.295 m，最大水深值为3.304 m；E6断面位于河道顺直段，水流较为平缓，水位较E5断面变化不大，最小水深值为3.359 m，最大水深值为3.303 m；E7断面位于弯道河段，河道发生较大的弯曲，惯性离心力的存在而使得自由水面的平衡状态遭到破坏，水流在弯道前的一段距离内，凸岸处自由水面略高于凹岸[15]，E7断面附近最小水深值为3.415 m，最大水深值为3.426 m；E8断面位于河道最下游，由于河道变窄，流量最大，因此该断面处的流速也最大，断面水位具有一定的变化，最小水深值为3.428 m，最大水深值为3.434 m。

图6 干流河道典型断面VOF分布

本研究根据河道起点、河道交叉点、河道顺直、河道急弯四个特点选取6个典型断面分析河道的水流速度分布，为河道的水生态分析提供指导，见图8。F1断面位于河道起点，河道断面较宽，宽度约29 m，由于该断面处于河道起点，无外源水补充，可以发现河道水流流速比较慢，且起点为盲端，不受外界水流扰动，水流流态非常稳定，表现为低流速下的旋流，水体平均流速为0.24 m/s；F2断面位于本研究河道与河道1的交叉口处，受河道补充水源的影响，该处水流表现为旋流、绕流、涡流、反射等紊流现象，表现为复杂的三维流动，水体平均流速为0.71 m/s；F3断面位于本研究河道与河道2的交叉口处，水流流态与F2断面类似，只是由于受河道2补水的影响，河道流量加大，水体流速变大，该处水体平均流速为0.83 m/s；F4断面位于本研究河道与河道3的交叉口处，该处河道边界条件非常复杂，一方面该处为两条河道的交汇处，另一方面该处为急弯断面，本研究河道的水流流向为自南向北，而河道3的水流流向为自北向南，所以两条河道的流向相反，由图8（F4）处的水流分析可以发现，河道水流整体还是按照本研究河道的流向流动，即自南向北，本研究河道有一部分水流流向河道3的下游，成为河道3的补充水源，而河道3的上游水流则主要汇入本研究河道，该处平均水流速度为0.84 m/s；F5断面位于本研究河道顺直段，河道中央水流较少受到外部边界条件影响，水流流态较为稳定，但是由于河道宽度不规则变化，在河岸处河水表现为较明显的紊流现象，并且受河岸阻力影响，水流流速较慢，该处平均水流速度为0.69 m/s；F6断面位于本研究河道急弯处，水流受河道转向的影响，在河岸处水体表现为明显的紊流现象，而中央水流则随河道转弯，该处平均水流速度为0.76 m/s。

图7 干流河道典型断面水深分布（单位：m）

图8 干流河道典型断面的流速分布（单位：m/s）

5 结 论

流水不腐，是减缓甚至基本消除河流黑臭的关键因素，因此对河道水动力特性进行数值模拟研究对消除河流黑臭具有重要意义。针对目前河流水动力数值模拟的研究多集中于一、二维，忽略了河道水流在竖直方向上的运动，无法精确模拟河流自由水面形状的复杂性和整体位置的不断变化的问题，本文采用CFD技术，考虑水流自由水面以及竖直方向的运动，建立了基于VOF法的河流水力三维紊流数学模型，并采用S.Soare Frazao与Y.Zech在急转弯河道模型中进行的试验结果对本模型进行了验证。以阜阳市某典型黑臭水体河道为研究对象，运用无结构贴体网格划分技术和局部加密网格生成技术建立了河道的计算网格模型，对河段水动力特性进行了三维数值模拟研究。结果表明：在河道受其他河流补水后，水深逐渐增大，河道最小水深为2.79 m，最大水深3.43 m；在河道起点盲段，水流为缓慢的旋流；在河道交叉以及急弯处，水流表现为旋流、绕流、涡流、反射等复杂的三维紊流现象，并且水流速度较顺直河段大，河道最小流速为0.24 m/s，最大流速为0.84 m/s。因此通过本文建立的三维紊流数学模型，实现了对自由水面、水深、流速的精确模拟，为河流水质的模拟提供基础，为河道水生态治理与恢复提供科学的决策支持。