堤防决口水流水力学特性的三维数值模拟

李火坤，邓冰梅，曾秀娟，张利荣，李宜忠

(1.南昌大学建筑工程学院，江西南昌 330031； 2.中国人民武装警察部队水电第二总队，江西南昌 330001)



堤防决口水流水力学特性的三维数值模拟

李火坤1，邓冰梅1，曾秀娟1，张利荣2，李宜忠2

(1.南昌大学建筑工程学院，江西南昌 330031； 2.中国人民武装警察部队水电第二总队，江西南昌 330001)

堤防决口水流的水力学特性对于制定决口封堵方案至关重要。基于FLOW-3D建立了模拟堤防决口水流的三维数值模型，首先对局部瞬间溃决水流水位场分布进行了三维数值模拟，模拟计算结果与C.Biscarini等代表性文献介绍的研究结果吻合较好，印证了数值模型的可靠性和准确性。以某流域河流堤防为例，建立了不同决口口门形状(U形、V形、T1形和T2形)的决口水流三维数值模型，计算分析了4种口门形状下决口附近的水位场及流速场分布规律。结果表明：进入决口后水流发生急速的侧向收缩和纵向跌落，决口处水深和流速在垂直水流方向断面上呈两侧小中间大的近似拱形分布，决口进口附近水流最大流速出现在决口左侧堤脚。数值模型和模拟方法可为制定堤防决口的最佳堵口方案提供技术参考。

堤防决口； 口门形状； 三维数值模拟； 水力学特性； FLOW-3D

堤防是指在江、河、湖、海沿岸或水库区、分蓄洪区周边修建的防止洪水漫溢或风暴袭击的挡水建筑物，是防洪体系中的重要组成部分。人类在发展过程中为兴水利、除水害而修建的堤防工程给社会带来巨大经济效益的同时存在着溃决风险，且在众多的洪涝灾害中，堤防决口是抢险最艰难、损失最严重、影响面最大的[1]。堤防溃决将直接威胁人民生命财产的安全，影响社会稳定。1953年荷兰风暴潮致900余处堤防溃决，1 835人死亡，直接经济损失达荷兰当年GDP的14%[2]；1998年长江流域中下游地区发生数十起严重的堤防溃决，大片田地、乡镇被淹[3-4]，造成巨大经济损失。堤防决口按照引发堤防决口的主要动力可分为水力决口型和非水力决口型两大类。影响堤防决口发展的因素主要可归纳为3个方面[5]：一是决口处的水动力因子，如决口处堤内外水位差、流速、负波加速度等；二是决口处的边界条件，如决口部位、大小、入流角、堤防土质结构特性以及含水率等；三是决口持续时间。

早期的堤坝溃决研究是以分析溃坝水流运动的理论解为主，1871年Saint-Venant 提出Saint-Venant 方程组后，J.J.Stoker[6]结合溃坝连续波和不连续波研究得到了溃坝波的Stoker 解；林秉南等[7]在假定河槽具有抛物线形断面且在平底、无阻力的简化条件下得出了有限长棱柱体水库的溃坝波对称解；谢任之[8]概括了现行不同的方法，较系统地提出了对于水库溃坝流量计算所归纳出的统一公式，给出了相关的表格；进入20世纪80年代后，数值模拟方面开始逐步广泛应用，R.Garcia等[9]将数学模型的数值解与基于Saint-Venant 方程组的一维无摩擦溃坝解析解进行了对比验证；R.J.Fennema[10]等采用不连续的边界条件模拟二维渐变流；C.Biscarini[11]等提出数值模拟溃坝水流自由表面流的浅水解决方案；魏文礼等[12]利用 Mac Cormack 预测-校正技术的隐式数值格式求解二维浅水方程，在此基础上建立模拟大坝溃决瞬间的洪水演进过程数学模型，预测矩形河道下游有多个障碍物时大坝局部溃决瞬间的洪水演进过程；冶运涛等[13]基于修正控制方程的复杂边界溃坝水流数值模拟，较好模拟复杂边界处水流特性；程永光[14]归纳总结了计算流体动力学最新理论知识，通过计算机和数值方法求解流体力学的控制方程，对流体力学问题进行模拟和分析。

在堤坝溃决后若能以决口现场的水力学特性参数为依据迅速有效地制定抢险措施进行封堵，将会大大减少受灾区人民的生命财产损失，本文以某流域河流堤防为例，基于FLOW-3D建立了不同决口口门形状(U形、V形、T1形和T2形)的堤坊决口水流三维数值模型，计算并分析了4种口门形状下决口附近的水位场及流速场分布，所建的数值模型和模拟方法可为制定最佳堵口方案提供技术参考。

1 基本理论

FLOW-3D为专业商业软件，采用独创的FAVOR法及真实的Tru-VOF法，具有离散完整的N-S方程，能很好地模拟真实世界中自由液面的流动现象并准确计算其流场。本次数值模拟选择非恒定、VOF及RNG湍流模型[15-16]。

1.1 基本方程

基本方程包括连续性方程和动量方程。

连续性方程：∂(uAx)/∂x+∂(vAy)/∂y+∂(wAz)/∂z=0

(1)

动量方程：

(2)

(3)

(4)

式中：u，v，w分别为x，y，z方向上的流速分量；Ax，Ay，Az为3个方向上可流动的面积分数；Gx，Gy，Gz为3个方向上的重力和非惯性力加速度；fx，fy，fz为3个方向上的黏滞力加速度；VF为可流动的体积分数；ρ为流体密度；p为作用在流体微元上的压强；t为时间。

1.2 RNG k-ε模型

RNG k-ε模型与标准k-ε模型相比，小尺度运动能从控制方程中去除，可以更好地处理流线弯曲较大及高应变率的流动。

(5)

(6)

1.3 自由表面处理技术

FLOW-3D软件对自由表面追踪采用VOF法。其原理[16]是模型中互不相融的两种或多种流体的任何相位，都由与其相对应的相关变量来描述，计算域整体的流动由所有相位的流动情况相加构成。任一单元的变量和特性取决于每一相的面积和体积分量值及其运动情况。

VOF运动学方程：

(7)

流体体积函数F=F(x,y,z,t)代表计算区域内流体的体积占计算区域的相对比例。F=1表示该单元完全被流体充满；F=0表示该单元完全被气体充满；0

2 决口水流数值模型算例模拟验证

图1 模型平面(单位：m)Fig.1 Plan of model (unit:m)

为模拟验证溃坝波的洪水演进过程，选取文献[9-10]中的局部瞬间溃决水位场计算算例进行验证。数值模型尺寸设为200 m×200 m×10 m(图1)，划分网格后的单元尺寸为2 m×2 m×1 m，总共1.0×105个单元。模拟的初始状态为坝上游水深10 m，下游水深5 m，坝中间缺口处不设置水体，底坡平顺且阻尼忽略不计。边界条件：设重力加速度为-9.8 m/s2，模型顶部设1个标准大气压，底部设为无滑移固壁边界，左右两侧皆设为固壁边界。求解设置：求解时间设为10 s，时间步长设为0.02 s。

在模拟时间为7.2 s时，溃坝波波前到达了河道左岸且在河道下游中心发展完全，因此，国内外学者都对其模型在7.2 s 时的水流计算结果进行分析比较，本文也取7.2 s时的水流计算结果(水深)与其他文献相比较。本算例所提取水深计算结果的截面A-A和B-B及点P1，P2的具体位置见图1[11]。将模拟计算结果与文献[11]中的二维和三维模拟及R.J.Fennema等[10]的散点结果进行对比(见图2)。

图2 水深计算结果与文献[11]计算结果对比Fig.2 Comparison between calculated results of depth given by paper and reference [11]

从图2可见，计算的水深数值大小及其变化规律与文献[11]中的二维及三维模拟结果基本吻合。所以，所建数值模型的计算结果与其他学者对该算例的计算结果吻合较好，验证了所建数值模型的可靠性以及对决口水流进行三维模拟的可行性。

3 堤防决口三维水流数值模型

堤防决口一般与天然河床水流方向垂直或斜交，决口呈突发性，除人工炸堤外，堤防决口一般都不是瞬间溃决，决口都有不同的口门扩宽过程，扩宽至一定程度后会趋于稳定。口门的扩宽受多种因素影响，扩宽过程十分复杂，为简化计算，对口门扩宽稳定后的堤防决口进行数值模拟计算，即采用瞬间溃决模式对不同口门形状下的决口水流进行模拟并分析其结果。以某实际流域河流堤防进行决口水流模拟建模，然后计算分析不同口门形状下的堤防决口水流的水力学特性。

图3 河道及堤防横断面(单位：m)Fig.3 Cross-sections of river channel and dike (unit:m)

模型河道内的模拟区域为长700 m，宽100 m(河底)的长直河道，模型高度10 m，河道坡降为0.3‰，河道糙率为0.025，决口高度8 m，距离下游河道200 m，河道左岸堤防高10 m，堤顶宽4 m，两侧边坡分别为1∶1.05和1∶0.85，河道及堤防典型断面如图3。此外，计算模型还包括决口处50 m内的平地段，综合考虑模型的尺寸和精度，将模型划分为3个网格块，第1个网格块为河道，网格划分为2 m×2 m×1 m；决口及决口出口50 m为另外2个网格块，对这2个网格块进行加密，将网格设为1 m×1 m×1 m。边界条件：物理边界条件设重力加速度为-9.8 m/s2，上游边界设流量为800 m3/s，水位为6 m，下游边界设为自由出流；模型顶部设1个标准大气压，模型底部和模型左右两侧皆设定为固壁边界。初始条件：给定河道初始水位为5 m。求解设置：求解时间设为600 s，时间步长设为2 s。

4 不同决口口门形状水流水力学特性的数值模拟

4.1 不同决口口门形状特征

实际决口的口门形状往往具有任意性，为简化计算，选择U形、V型、T1形和T2形这4种典型的决口口门形状进行数值模拟。不同口门形状决口模型的顶部长20 m，宽30 m，高8 m，具体形状及尺寸如图4所示。经过计算得到不同口门形状下各决口形状水流在t=523 s后达到稳定，选取t=523 s时刻的计算数据进行分析。

图4 不同形状决口口门(单位：m )Fig.4 Different shapes of dike breachs (unit：m)

图5 决口处平面截面示意(单位：m )Fig.5 Dike breaking sections (unit：m)

4.2 水位场分布规律

选择决口水流C-C，D-D，E-E，F-F等截面计算结果进行分析，各截面位置如图5所示。计算结果表明，总体上不同口门形状的决口处水流的水深分布规律基本相似，水流在进入决口后发生急速的侧向收缩，在沿决口水流方向(y方向)，水流距离堤脚处越远，水位越低；在决口宽度方向(x方向)，两侧堤脚处特别是靠近河道上游来水侧水深较低，而口门中部位置相对较高，水面线在垂直于决口中轴线方向上呈两边低、中间高的拱形，水流在决口两侧堤脚处产生水跌，形成卷吸；各不同决口口门形状的水深计算结果见图6和7。

图6 不同形状决口水深计算结果三维图Fig.6 3D depth hydrographs for dike breakdown with different shapes

图7 不同断面水深分布Fig.7 Depth hydrographs along different sections

图7(a)为决口C-C断面水面线沿水流方向的分布，图7(b)～(d)为决口D-D,E-E,F-F等断面水面线沿决口宽度方向的分布。不同口门形状的决口水流水深沿水流方向变化规律基本一致，均沿水流方向水面迅速降低。在决口进口D-D断面，4种口门形状的断面水深整体较高，水深总体上呈两边低中间高的分布特征，且决口右侧水深较左侧水深高；各口门形状在该断面水深最大值由大到小分别为T1形(x=18.5 m，水深3.87 m)、V形(x=14.5 m，水深3.69 m)、U形(x=19.5 m，水深3.36 m)、T2形(x=23.5 m，水深3.33 m)。决口中部的E-E断面和决口出口F-F断面的水面线呈中间较高的近似拱形分布，由于V形决口口门宽度最小，故该形状决口的水面线最窄，最高水位出现在接近决口中轴线附近；其他形状决口在该两个断面的最高水位多数出现在决口右侧位置，且随决口断面面积增大其最高水位越靠近决口右边界处。

4.3 流速场分布规律

各口门形状的决口流速计算结果见图8。

图8 不同形状决口流速场三维图和流速矢量图(单位：m·s-1)Fig.8 3D velocity diagrams and velocity vectors for velocity field having different shapes (unit:m·s-1)

由图8可见，这4种口门形状的决口流速分布特征基本相似，在决口处沿水流方向，距堤脚越远，流速越大；在决口进口横断面，最大流速出现在决口左侧，由于部分上游水流发生转向，决口出口处又无水流的顶托作用，因此流速剧增；在决口出口断面，水流向四周平地扩散，水面跌落较大，产生较大流速。

决口处断面流速分布见图9。

图9 不同断面流速分布Fig.9 Velocity distributions along different sections

由图9可见，C-C断面上，流速沿水流方向增大；在决口宽度方向，4种口门决口处的流速分布基本相似，均呈不对称拱形分布。在进口D-D断面，最大流速均出现在决口中轴线左侧，U形决口在该断面流速最大，最大流速4.40 m/s，出现位置为x=6.5 m，V形决口在该断面流速最大值中最小，最大流速为2.45 m/s，出现位置为x=12.5 m。决口E-E和F-F断面分别靠近下游和出口位置，过流断面有所减小，流速较进口断面大，流速沿宽度方向的分布呈中间大两侧小的拱形分布；对于E-E断面，除V形决口外，其他口门形状决口最大流速出现在决口左侧位置；对于决口出口的F-F断面，流速在3个断面中最大，各口门形状流速最大值均在决口中轴线右侧。

5 结 语

计算分析了在上游流量和水头及口门宽度一定的条件下不同口门形状的决口处水位场和流速场分布特征。总体上不同口门形状的决口处水深、流速分布基本相似，数值上有所差别。各断面中水深最大值为T1形决口在x=18.5 m处，水深为3.87 m；决口进口断面最大流速均出现在决口中轴线左侧，最大值为U形决口在x=6.5 m处，流速为4.40 m/s。对决口封堵而言，决口处的流速分布对科学制定决口封堵方案及选择封堵料粒径尤为重要，从上述各决口口门形状的流速分布特征来看，决口口门左侧流速较大，极易造成该侧堤防堤脚的冲刷从而扩大决口宽度，因此在封堵决口时，应及时对该侧决口堤脚进行裹头处理；在选择封堵料粒径时，左侧口门处的封堵粒径要比右侧口门封堵粒径大才可保证该侧封堵料的抛投稳定。本文所建的决口水流水力学特性模拟的三维数值模型和模拟方法可为堤防决口的封堵提供技术支持和参考。

[1]张永忠.抗洪抢险技术[M].北京：军事科学出版社,1999.(ZHANG Yong-zhong.Flood rescue technology[M].Beijing：Military Science Publishing House,1999.(in Chinese))

[2]ZHU Y H,VISSER P J,VRIJLING J K,et al.Experimental investigation on breaching of embankments[J].Sci China Tech Sci,2011,54(1)：148- 155.

[3]张我华,方仲将,任延鸿.九江大堤溃坝渗流的可靠性分析[J].浙江大学学报(工学版),2005,39(7)：976- 982.(ZHANG Wo-hua,FANG Zhong-jiang,REN Yan-hong.Reliability analysis for seepage breach of Jiujiang dike on Changjiang River[J].Journal of Zhejiang University(Engineering Science),2005,39(7)：976- 982.(in Chinese))

[4]杨光煦.河堤决口水力特性及堵口技术[J].湖北水力发电，1999,35(4)：9- 16.(YANG Guang-xu.Hydraulic characteristics of breached river dike and plugging technology[J].Hubei Water Power,1999,35(4)：9- 16.(in Chinese))

[5]曲红玲.河道溃堤与溃堤波的一、二维耦合计算数值模拟[D].南京：河海大学,2006.(QU Hong-ling.Numercial simulation for coupling one-dimensional river and two-dimensional dike-break waves[D].Nanjing：Hohai University,2006.(in Chinese))

[6]STOKER J J.Water wave[M].New York：Interscience Publishers,1957：334- 341.

[7]林秉南,龚振瀛,王连祥.突泄坝址过程线简化分析[J].清华大学学报,1980,20(1)：17- 31.(LIN Bin-nan,GONG Zhen-ying,WANG Lian-xiang.Dam-site hydrographs due to sudden release[J].Journal of Tsinghua University,1980,20(1)：17- 31.(in Chinese))

[8]谢任之.溃坝水力学[M].济南：山东科学技术出版社,1993.(XIE Ren-zhi.Dam break hydraulics[M].Jinan：Shandong Science and Technology Press,1993.(in Chinese))

[9]GARCIA R,KAHAWITA R A.Numerical solution of the St.Venant equations with the MacCormack finite-differences scheme[J].International Journal for Numerical Methods in Fluids,1986,6(5)：259- 274.

[10]FENNEMA R J,HANIF CHAUDHRY M.Explicit methods for two-dimensional transient free-surface flows[J].Journal of Hydraulic Research,1990,116(8)：1013- 1034.

[11]BISCARINI C,FRANCESCO S D,MANCIOLA P.CFD modelling approach for dam break flow studies[J].Hydrology & Earth System Sciences,2010,14(4)：705- 718.

[12]魏文礼,沈永明,孙广才,等.二维溃坝洪水波演进的数值模拟[J].水利学报,2003(9)：43- 47.(WEI Wen-li,SHEN Yong-ming,SUN Guang-cai,et al.Numerical simulation of 2D dam-break flood wave[J].Journal of Hydraulic Engineering,2003(9)：43- 47.(in Chinese))

[13]冶运涛,梁犁丽,张光辉,等.基于修正控制方程的复杂边界溃坝水流数值模拟[J].水力发电学报,2014,33(5)：99- 107.(YE Yun-tao,LIANG Li-li,ZHANG Guang-hui,et al.Numerical simulation of dam-break water flow with complex boundary based on governing equations modification[J].Journal of Hydroelectric Engineering,2014,33(5)：99- 107.(in Chinese))

[14]程永光.计算流体动力学[M].北京：中国水利水电出版社,2014.(CHENG Yong-guang.Computational fluid dynamics[M].Beijing：China Water Power Press,2014.(in Chinese))

[15]张婷.波浪的三维数值模拟及其应用[D].天津：天津大学,2009.(ZHANG Ting.Three-dimensional numerical simulation of waves and its application[D].Tianjin：Tianjin University,2009.(in Chinese))

[16]张健,方杰,范波芹.VOF方法理论与应用综述[J].水利水电科技进展,2005,25(2)：67- 70.(ZHANG Jian,FANG Jie,FAN Bo-qin.Advances in research of VOF method[J].Advances in Science and Technology of Water Resources,2005,25(2)：67- 70.(in Chinese))

3D numerical simulation of hydraulics characteristics of levee breach flow

LI Huo-kun1，DENG Bing-mei1，ZENG Xiu-juan1，ZHANG Li-rong2，LI Yi-zhong2

(1.SchoolofCivilEngineeringandArchitecture,NanchangUniversity,Nanchang330031,China; 2.TheNo.2GeneralTeamoftheChineseArmedPoliceHydropowerTroops,Nanchang330001,China)

The hydraulic characteristics of flow caused by the levee breach are very important for working out the schemes for closure of breach when there is a sudden levee breach due to breaking.A 3D numerical model is established to simulate the sudden levee breach flow based on the FLOW-3D software,by which,first,the 3D numerical simulation of distribution of the water level field induced by local instantaneous breaking flow has been carried out.And it is found that the calculated results based on FLOW-3D software have a good agreement with the representative research results given by C.Biscarini,which has verified the reliability and accuracy of the numerical model in this paper.Taking the levee located on a river basin as an example,a 3D numerical model is set up for the levee breakdown water flow outflowing from different levee breach ( including U-shaped,V-shaped,T1 and T2-shaped breaches).At the same time,and the distribution characteristics of the water level and velocity field near the breach in the case of these four breach shapes were calculated and analyzed.The simulated and analyzed results show that there was a rapidly lateral contraction and longitudinal fall of the water flow after inflowing into the breach,and the water depth and velocity lines at the breach appear to be approximate arched profile,which being small on both sides and large in the middle along the vertical water flow cross-section direction.And the maximum water flow velocities close to the levee breach occurred at the left toe.The numerical model and the simulation methods based on FLOW-3D software presented by this paper can provide a technical reference for working out the best schemes for the emergency closure of breach.

levee breach; levee breach shapes; 3D numerical simulation; hydraulic characteristics; FLOW-3D software

10.16198/j.cnki.1009-640X.2016.06.005

李火坤，邓冰梅，曾秀娟，等.堤防决口水流水力学特性的三维数值模拟[J].水利水运工程学报,2016(6):29-36.(LI Huo-kun,DENG Bing-mei,ZENG Xiu-juan,et al.3D numerical simulation of hydraulics characteristics of levee breach flow[J].Hydro-Science and Engineering,2016(6):29-36.)

2015-11-05

国家自然科学基金资助项目(51269019，51469015) ；广东省水利科技创新基金资助项目(2014-08)

李火坤(1981—)，男，湖南长沙人，博士，副教授，主要从事水工水力学及泄流结构动力检测方面的研究工作。E-mail:lihuokun@126.com

TV871； TV122+.4

A

1009-640X(2016)06-0029-08