基于Flow-3D的堤防堵口水力条件数值模拟

王 兵，张凤德，贺昌海

(1. 武汉大学水资源与水电工程科学国家重点实验室，武汉 430072；2. 黑龙江省三江工程建设管理局，哈尔滨 150081)

堤坝溃决水流一般兼有急流和缓流、涌波和稀疏波等流态，决口流量过程常是峰高量大，形状陡峭，急缓流转换，水面梯度较大，在口门区域常存在间断[1,2]。溃堤问题与溃坝问题相比最大的区别在于河道的侧向水流运动，与侧堰水流运动非常相似[3]。目前针对河道堤防溃决水流的相关数值模拟研究，主要为二维平面数值模拟[4-6]。

为了尽量减少堤防溃决后造成的损失，结合实际地形地质条件[7-9]和相应口门区域水力条件及变化规律，制定科学合理的堵口抢险预案在当前时代特征下显得尤为重要。在传统物理模型试验存在周期长、投资大等缺点的背景下，近年来Flow-3D等CFD软件得到了越来越广泛的运用，包括溢流坝泄流[10]、消能池水跃[11]、泥沙输运[12]、长廊式气垫调压室的水-气两相流[13]等三维数值模拟计算。众多前人的研究成果表明，三维水流数值模拟计算结果与物理模型试验结果是相对一致的[14]，可用于指导和解决实际工程问题。

本文利用Flow-3D软件，针对某河道砂性堤防溃决水流运动进行三维数值模拟，研究超标准洪水条件下口门区域水力条件，其成果可为研究极端洪水条件下该河道砂性堤防溃口快速封堵技术提供科学依据。

1 研究河道概况

某河道两岸地形地貌为低漫滩阶地，滩地高程多在41～44 m之间，而岛面高程多为41～42 m。河床与河岸的土质以细砂及壤土为主，抗冲刷能力差。平槽水位时河道水深达10 m以上，水面宽阔，水深流急。该河道堤防防洪标准为50 a一遇，堤防形式以砂基砂堤或砂基土堤为主，部分堤段受地形地质等因素的影响为险工弱段。

研究选取某个典型堤防险段及其对应的河道实际地形，对百年一遇洪水条件下的堤防溃口水流运动进行数值模拟。该典型险段设计断面如图1，堤顶超高2.2 m，堤顶宽8 m，两侧边坡均为1∶4.0。

图1 典型险段设计断面(单位：cm)Fig.1 Design section of typical risk dike

2 数值计算

2.1 控制方程

本文研究的具有实际复杂地形的堤防溃口水流问题，因地势起伏多变及河道的侧向水流运动特性，溃口区域水流流态复杂、水面梯度大，因此采用RNGk-ε紊流模型。控制方程如下：

连续方程：

(1)

动量方程：

(2)

紊动能k方程：

(3)

紊动能耗散率ε方程：

(4)

式中：ui是流速张量分量；Ai、gi、fi分别代表三维直角坐标方向上可流动的面积分数、重力加速度和黏滞力；VF是可流动的体积分数；ρ是流体密度；p是作用在流体微元上的压力；k、ε分别是紊动能和紊动能耗散率；σk、σε分别是紊动能和耗散率对应的Prandtl数；μ、μt分别是水的动力黏滞系数和紊动黏滞系数；Gk是紊动能k的产生项；C*ε1、Cε2是经验常数。

采用交错矩形网格的FAVOR法离散控制方程，采用TruVOF方法处理流体自由表面，采用极小残差(GMRES)算法求解压力。

2.2 模型及网格划分

在原有地形资料的基础上，利用CATIA软件进行三维河道地形建模，并创建堤防建筑物。将生成的模型实体导入Flow-3D中进行数值计算。为了保证计算的效率、稳定性和准确性，采用均匀的立方体网格，纵横比为1；并考虑深度方向的计算精度要求相对较高，设置网格尺寸X∶Y∶Z=2∶2∶1。为更精确的模拟溃口区域的水力条件，在溃口区域采用嵌套网格。网格块①覆盖河道范围对应堤线桩号从1+000 m至3+150 m处(原型)，立方体网格尺寸为4 m×4 m×2 m；嵌套网格块②覆盖范围为溃口口门堤线前后200 m、上下堤头各退后50 m所围成的区域，溃口中心对应堤线桩号2+100 m处，立方体网格尺寸为2 m×2 m×1 m。因此网格总数随口门宽度变化而稍有变化，例如口门宽60 m时，网格总数约379万，有效网格数约为229万。模型及网格划分如图2所示，(a)中网格单元扩大六倍以方便显示。

图2 模型及网格划分示意图Fig.2 Schematic diagram of model and grid division

2.3 边界条件及初始条件

边界条件(Mesh①)：上游Ymin设为流量边界，根据所建模型地形断面数据和百年一遇洪水位46.17 m换算而得；下游Ymax设为压力边界，用下游水位控制，下游水位通过试算确定；左岸Xmin边界与流体不接触，不影响流态，采用缺省设定为对称边界；溃口堤后Xmax设为压力边界，根据与百年一遇洪水位(溃口处河道断面)的水位落差ΔZ1控制；模型底部Zmin设为固壁边界；模型顶部Zmax设为压力边界，P=0(一个大气压)。Mesh②除顶部与底部外，均取对称边界。

初始条件：以下游控制水位作为计算域内的初始水位。

流量监测：Baffles 在FLOW-3D中是没有厚度的孔隙孔板，用以控制或引导水流。实际工作中常用于测量通过某断面的水流流量与计算通过的颗粒数量，且完全不影响水流运动。因此设置Baffle1、Baffle2和Baffle3分别监测模型上游进口流量Q1、下游出口流量Q2和溃口处分流量QK，如图2所示。

2.4 计算工况

研究主要利用FLOW-3D分别对不同口门宽度及堤防内侧不同水位条件作用下的溃口水流运动进行数值模拟，得到不同工况条件下口门区域的水位、水位场分布、流速及流速场分布等水力学特性，进而分析其对口门区域水力条件的影响程度。计算工况如表1。

表1 数值模拟计算工况Tab.1 Calculation conditions of numerical simulation

3 数值计算结果分析

3.1 不同口门宽度下流速场

图3为沿堤防轴线口门断面的沿水深平均流速分布图，在不同口门宽度下其流速分布形态基本一致。口门断面上，靠近两堤头流速急剧减小，中间段流速分布较堤外断面更为均匀，但下堤头前出现局部流速峰值，而数值随口门束窄先增大后减小。图3(b)与图3(a)相比可知，随着ΔZ1减小，也即堤后水位上涨，口门宽度相同条件下断面流速明显减小。

图3 口门断面沿水深平均流速Fig.3 Depth-averaged velocity of the entrance section

堤外流速峰值分布于口门偏上游区域，随着口门的束窄流速变化坡度增大。堤内流速集中于口门中间区域，横向扩散现象明显。例如口门宽60 m时口门区域流速场，如图4。

图4 口门区域流速场Fig.4 Velocity field of the entrance area

3.2 不同口门宽度下水位场

图5为沿口门中轴线断面水面线分布图，溃口水流水面梯度随口门宽度的减小而增大。堤线前水位随口门宽度减小而升高，水流经过口门过程中出现垂向收缩的现象。

图5 口门中轴线水面线Fig.5 Water surface profile of the central axis of the entrance

在上、下堤头前，堤角附近的卷吸现象明显，对口门区域横向水位变幅影响较大。例如口门宽60 m时口门区域水位场，如图6。

图6 口门区域水位场Fig.6 Free surface elevation of the entrance area

3.3 局部大流速分布特征

根据数模计算结果，堵口进占过程中口门区域存在局部大流速分布。60 m龙口宽度时口门断面流速分布如图7，下堤头附近流速明显高于上堤头，工况7条件下最大流速达8.27 m/s，工况8条件下最大流速达6.74 m/s，位置向下堤头前靠拢；同时口门底部流速较大。因此合龙前必须分别针对下堤头和口门底部的冲刷情况采取相应的裹护措施。

图7 口门断面流速分布图Fig.7 Velocity distribution diagram of breach section

3.4 溃口分流量和口门断面平均流速

根据Baffles监测结果，堤防溃口口门宽度与堤后边界条件对溃口分流量的影响十分显著。各工况条件下溃口分流量与口门断面平均流速如表2。

4 结 语

研究选取某河道堤防典型险段实际地形和现有防洪设计标准(P=2%)的堤防建筑物，采用CATIA软件作为三维建模平台，并基于FLOW-3D软件对百年一遇洪水条件下、不同口门宽度及堤防内侧不同水位条件作用下的溃口水流进行数值模拟，从而得到14种组合工况条件下口门区域的水位、水位场分布、流速及流速场分布等堤防堵口水力学特性。计算结果表明，溃口内外水位差和口门宽度均对溃口分流量、口门断面平均流速有显著的影响，应当据此选择合适的堵口时机；随着口门束窄，溃口前水位逐渐抬高，口门断面流速先增后减，堵口进占过程中因根据流速变换抛投材料及粒径；合龙过程中应适时监测口门底部和下堤头前的冲刷现象，并采取相应的防护措施；为保证堵口截流顺利实施，应当具备足够的抛投强度等。

表2 溃口分流量和口门断面平均流速Tab.2 Flow rate and average velocity of breach section

注：堤身内外实际落差ΔZ2以口门中轴线上距离原堤轴线100 m的两点处水位差计。

□

[1] 谭维炎．计算浅水动力学[M]．北京：清华大学出版社，1998．

[2] 张大伟．堤坝溃决水流数学模型及其应用研究[D]．北京：清华大学，2008.

[3] 张修忠．堤防溃决过程的数值模拟[D]．北京：清华大学，2003.

[4] 苑希民,田福昌,王丽娜．漫溃堤洪水联算全二维水动力模型及应用[J].水科学进展，2015，26(1)：83-90.

[5] 贺治国.考虑水沙相互作用的河堤溃决二维数值模型[J].水动力学研究与进展A辑，2010，25(2)：147-154.

[6] 陈 珺，张小峰，谈广鸣，等.考虑溃口展宽的溃堤水流泥沙数值模拟[J].水动力学研究与进展A辑，2007，22(5)：647-653.

[7] 杨佩国，杨勤业，吴绍洪,等．基于数值模拟的黄河下游不同情景溃堤洪水特性[J]．地理研究，2007，26(2)：328-336．

[8] 张大伟，张 超，王兴奎．具有实际地形的溃堤水流数值模拟[J]．清华大学学报(自然科学版)，2007，47(12)：2 127-2 130．

[9] 徐国宾，孟庆林，苑希民．含沙溃堤洪水数值模拟分析[J]．天津大学学报(自然科学与工程技术版)，2016，49(10)：1 008-1 015．

[10] 戚 蓝，陈 辉，费文才，等．基于精确河道地形的溢流坝泄流三维数值模拟[J]．水力发电学报，2016，35(10)：12-20．

[11] 王月华，包中进，王 斌．基于Flow-3D软件的消能池三维水流数值模拟[J]．武汉大学学报(工学版)，2012，45(4)：454-457．

[12] 刘成林，陈宇豪．基于Flow-3D的水平射流冲刷泥沙数值模拟[J]．人民长江，2016，47(6)：87-91．

[13] 夏林生，程永光，周大庆．3-D simulation of transient flow patterns in a corridor-shaped air-cushion surge chamber based on computational fluid dynamics[J].Journal of Hydrodynamics，2013，25(2)：249-257.

[14] 贺昌海，陈 辉，刘 全．基于CATIA建模的导流工程三维数值模拟研究[J]．天津大学学报(自然科学与工程技术版)，2016，49(4)：422-428．