涌潮作用下丁坝水力特性的数值模拟研究

何 昆，杨元平，陈 刚，钱盛杰，张芝永

（1.浙江省水利河口研究院（浙江省海洋规划设计研究院）， 浙江 杭州 310017；2.海宁市盐仓江堤管理所， 浙江 海宁 314400）

1 问题的提出

丁坝作为一种护岸建筑物与河道整治建筑物，在水利工程中应用广泛。涌潮是一种发生于喇叭形河口或海湾地区的特殊潮汐现象，具有潮差大、水动力强的特点，易对丁坝造成冲刷损坏。因此，研究涌潮作用下丁坝的水力特性，对丁坝的布置优化及安全运行具有重要意义。

对丁坝水力特性的研究一直是河流动力学科的热门研究方向之一，研究方法主要有物理模型试验和数值模拟等。1953 年Ahmed[1]首次对丁坝绕流的水面线变化情况进行测量；彭静[2]等基于颜料示踪和油膜技术对丁坝群附近的流场分布进行可视化记录研究；Jennifer[3]等采用激光多普勒测速仪（LDV）对丁坝附近的水位和流速分布进行观测，并在此基础上研究丁坝水流的三维紊动特性；凌建明[4]等运用二维水深平均k-ε 紊流模型，采用有限元的方法对控制方程进行离散，通过计算得到不同工况下绕流丁坝附近的流场分布；刘超[5]等采用有限体积法离散求解二维浅水方程，辅以RNG k-ε湍流模型，模拟丁坝坝头附近的湍流流场分布；胡子俊[6]等采用RANS 方法对丁坝绕流进行模拟，将不同湍流模型与试验观测结果进行对比后得出RNG k-ε 湍流模型的精度较高，通过数值计算得到的流场变化形态与试验观测和客观流场运动规律吻合的结论。

但目前对丁坝的三维数值模拟或物理模型模拟，大部分针对恒定水流形式，对涌潮作用下丁坝水力特性的研究较少。本文通过建立三维数值水槽模型，对涌潮作用下丁坝的水力特性规律进行探索，结果可为涌潮作用地区的丁坝设计提供参考。

2 理论模型简介

2.1 控制方程

模拟的水体运动为不可压缩黏性流体在笛卡尔坐标系中的运动，其运动的控制方程为连续性方程及动量守恒方程（Navier-Stokes 方程）。其中，由于采用VOF 算法进行水体自由表面追踪，在连续性方程和动量方程中加入体积和面积分数参数，具体表达式如下：

连续性方程：

动量方程：

式中：ρ为流体密度，kg/m3；p为压强，Pa；VF为流体的体积分数；Ax，Ay，Az是x，y，z3 个方向流体的面积分数；u，v，w为x，y，z3 个方向的速度分量，m/s；Gx，Gy，Gz为x，y，z3 个方向的重力加速度，m/s2；fx，fy，fz为在x，y，z3 个方向的流体黏滞力加速度，m/s2。

2.2 紊流模型

流体与丁坝作用后通常会产生变形破碎和湍流漩涡，为更好地模拟水体运动状态，本文采用RNG k-ε 模型封闭水流运动控制方程。其方程式如下所示：

式中：kT为紊动动能，J；PT为紊动动能的表征项；GT是浮力作用下的紊动动能产生项，在不可压缩流体中默认为0；εT为紊流能量消散率；CDIS为无量纲参数，DiffT和DiffS为能量消散项。

2.3 方程的离散与求解

采用有限差分法进行数值离散，将空间计算域离散成若干个结构化矩形网格单元求解方程。计算网格为交错网格，压强、流体体积分数、涡黏滞系数、湍流动能、耗散率等标量均定义在控制体的中心，速度等矢量则定义在垂直于矢量方向的控制体表面中心。

方程离散后，选用GMRES 迭代算法进行压力迭代求解。动量方程中的求解，本文选择显式解法中二阶精度格式进行离散，该格式能够更好地模拟水流的变形破碎。

2.4 模拟范围与网格划分

计算模型范围在X方向（水流方向）以丁坝为中心的长度取20L，Y方向取5L，Z方向取H1。X、Y方向采用非均匀矩形网格划分，并对丁坝附近网格进行局部加密计算。其中非加密部分X、Y方向网格尺度最大取0.2L， 加密部分X、Y方向网格尺度最大取0.1L（见图1）。Z方向网格为均匀网格，网格尺度取0.1H1。（L为坝长，H1为坝高）。

图1 XY 平面网格划分示意图

2.5 边界条件

数值模拟水池左侧即Xmin使用速度边界条件，通过给定边界流体的流速和水位来模拟生成恒定流以及涌潮；右侧Xmax使用出流边界；顶端Zmax使用压力边界条件，底部Zmin为墙边界条件；前侧Ymin为墙边界条件，后侧Ymax为对称边界条件。

3 模型验证

涌潮作用下的丁坝水力特性物理模型试验存在参数测量难度大且精度不高的问题，因此，本文采用恒定流作用下的丁坝水流结构进行数学模型验证。验证试验在浙江省河口海岸重点实验室的玻璃水槽中进行。水槽长50.0 m，宽1.2 m，高0.6 m。试验所用丁坝长度L=40 cm，于水槽中间段中心位置对称布置，试验水深d=20 cm。丁坝由厚度为10 mm的有机玻璃制成，其高度取22 cm，略高于水位以模拟非淹没式刚性坝。试验立面示意见图2。

图2 试验立面示意图

试验水流条件见表1。流速测量采用“小威龙”剖面流速仪，纵向流速测点断面布置依据距离结构物近则密、远则疏的原则，并考虑流场对称性。为提高效率，横向取1/2 结构布置测点，从中心点起沿水流方向（纵向）布置0～4 共5 个测流断面，断面间距均为10 cm。每个纵向断面的上游布置S1～S3 共3 个断面，下游布置0～9 共10 个断面，断面间距及具体布置见图3。

表 1 试验水流条件表

图3 流速测点平面布置图 单位：cm

各工况下，表层断面（z=18 cm）处试验测量结果与数值模拟结果的速度云图对比见图4。从图4 中可以看出，数值模拟与试验的流速分布情况及规律基本相同，由于丁坝的挑流作用及过水断面缩窄，坝头附近流速均显著增大。同时由于丁坝的掩护作用，下游丁坝正后方区域流速显著减小。由于试验测量断面数据插值与数值模拟数据插值相比较为稀疏，因此坝头附近的流速与实测流速有一定差异，坝下游侧的回流区也有一定误差。总体而言，本文建立的数学模型能够相对准确地模拟丁坝影响下的水力特性。

图4 表层断面速度云图对比图

4 水力特性分析

采用参考文献7 中的方法及数值生成涌潮，并在此基础上模拟涌潮作用下非淹没丁坝附近的水力特性情况。以潮前水深d=10 cm、潮头高度H=5 cm 及传播速度c=0.45 m/s 的涌潮为例分析说明。

4.1 水位分析

涌潮作用下的水位分布情况见图5。由图5 可知，涌潮在潮头刚刚到达丁坝时，由于丁坝的阻流效应，迎潮侧水位液面快速抬高，最大壅高可达2H。背潮侧水位未发生明显变化，仅接近丁坝的潮头线明显晚于其他位置处的潮头线，且越接近坝根处，潮头线的滞后时间越长。

图5 涌潮作用下水位分布图

涌潮潮头经过丁坝后，迎潮侧水位壅高明显减小，最大壅高仅0.8H。涌潮作用下，背潮侧的回流区范围内水位下降明显，最大下降幅度达H。

4.2 流速分析

涌潮作用下的平均流速分布情况见图6。由图6 可知，当涌潮潮头刚刚到达丁坝时，丁坝的束流效果较弱，仅在沿丁坝长度方向L范围内的水流产生绕流，流向发生变化。流速大小发生变化的范围基本在丁坝附近，坝头附近局部范围内流速增大，最大流速达1.7c。丁坝迎潮侧水流受丁坝壅阻，流速减缓。背潮侧流速分布基本以潮头线为界，潮头线至丁坝范围内流速为涌潮传播速度，潮头线至丁坝范围外流速接近于0。

图6 涌潮作用下平均流速分布图

在涌潮潮头经过丁坝后，丁坝的束流效果增强，沿丁坝长度方向4L范围内的水流流向发生变化，且在丁坝背潮侧产生回流漩涡。坝头背潮侧附近流速增大的范围相比刚到达时大幅增加，沿涌潮传播方向长度达L，垂直涌潮传播方向长度达1.2L，该范围内的最大流速略减小至1.5c，且丁坝坝根附近迎背侧均存在流速接近于0 的区域。

4.3 近底切应力分析

涌潮作用下的近底切应力分布情况见图7。由图7 可知，在涌潮潮头刚刚到达丁坝时，切应力与流速分布相似，背潮侧及靠近丁坝1.5L范围内的迎潮侧，切应力接近于0，其余位置处的切应力达0.4 Pa，仅坝头位置局部范围内的切应力较大，近底最大切应力达1.7 Pa。

图7 涌潮作用下近底切应力分布图

涌潮潮头经过丁坝后，在沿丁坝长度方向L范围内的切应力基本都接近于0。除坝头背潮侧位置沿涌潮传播方向长度达L，垂直涌潮传播方向长度达1.5L的范围内切应力相对较大，平均切应力达1.1 Pa，最大紊动切应力达1.9 Pa，其余位置处的切应力基本在0.4 Pa 以内。由此可见，坝头背潮侧位置切应力相对较大，泥沙易在此起动形成冲刷作用，而迎潮侧则由于切应力始终较小，水流携带的泥沙易在此沉降形成淤积作用。

5 结 语

在利用物理模型试验成果对数学模型进行验证的基础上，模拟研究了涌潮作用下非淹没丁坝的水位、流速、切应力分布情况，得到以下主要结论：

涌潮潮头刚到达丁坝时，多表现为波浪的特性，水位在丁坝迎潮侧有明显壅高，流速及切应力分布以潮头线为界，仅坝头局部范围内扰流数值较大；涌潮潮头经过丁坝后，多表现为水流的特性，丁坝迎潮侧水位壅高，背潮侧水位下降，流速及切应力在坝头背潮侧沿涌潮传播方向L，垂直涌潮传播方向长度1.5L的范围内较大，这也是涌潮作用下丁坝坝头冲刷需要重点防护的区域。

由于实际工程中丁坝的结构形式及布置方式各有不同，本次模拟的水力特性与实际有一定差异，对丁坝群等与工程实际更为符合的情况有待进一步研究。