李火坤，鄧冰梅，曾秀娟，張利榮，李宜忠

(1.南昌大學建筑工程學院，江西南昌 330031； 2.中國人民武裝警察部隊水電第二總隊，江西南昌 330001)

?

堤防決口水流水力學特性的三維數(shù)值模擬

李火坤1，鄧冰梅1，曾秀娟1，張利榮2，李宜忠2

(1.南昌大學建筑工程學院，江西南昌 330031； 2.中國人民武裝警察部隊水電第二總隊，江西南昌 330001)

堤防決口水流的水力學特性對于制定決口封堵方案至關重要?；贔LOW-3D建立了模擬堤防決口水流的三維數(shù)值模型，首先對局部瞬間潰決水流水位場分布進行了三維數(shù)值模擬，模擬計算結(jié)果與C.Biscarini等代表性文獻介紹的研究結(jié)果吻合較好，印證了數(shù)值模型的可靠性和準確性。以某流域河流堤防為例，建立了不同決口口門形狀(U形、V形、T1形和T2形)的決口水流三維數(shù)值模型，計算分析了4種口門形狀下決口附近的水位場及流速場分布規(guī)律。結(jié)果表明：進入決口后水流發(fā)生急速的側(cè)向收縮和縱向跌落，決口處水深和流速在垂直水流方向斷面上呈兩側(cè)小中間大的近似拱形分布，決口進口附近水流最大流速出現(xiàn)在決口左側(cè)堤腳。數(shù)值模型和模擬方法可為制定堤防決口的最佳堵口方案提供技術參考。

堤防決口； 口門形狀； 三維數(shù)值模擬； 水力學特性； FLOW-3D

堤防是指在江、河、湖、海沿岸或水庫區(qū)、分蓄洪區(qū)周邊修建的防止洪水漫溢或風暴襲擊的擋水建筑物，是防洪體系中的重要組成部分。人類在發(fā)展過程中為興水利、除水害而修建的堤防工程給社會帶來巨大經(jīng)濟效益的同時存在著潰決風險，且在眾多的洪澇災害中，堤防決口是搶險最艱難、損失最嚴重、影響面最大的[1]。堤防潰決將直接威脅人民生命財產(chǎn)的安全，影響社會穩(wěn)定。1953年荷蘭風暴潮致900余處堤防潰決，1 835人死亡，直接經(jīng)濟損失達荷蘭當年GDP的14%[2]；1998年長江流域中下游地區(qū)發(fā)生數(shù)十起嚴重的堤防潰決，大片田地、鄉(xiāng)鎮(zhèn)被淹[3-4]，造成巨大經(jīng)濟損失。堤防決口按照引發(fā)堤防決口的主要動力可分為水力決口型和非水力決口型兩大類。影響堤防決口發(fā)展的因素主要可歸納為3個方面[5]：一是決口處的水動力因子，如決口處堤內(nèi)外水位差、流速、負波加速度等；二是決口處的邊界條件，如決口部位、大小、入流角、堤防土質(zhì)結(jié)構(gòu)特性以及含水率等；三是決口持續(xù)時間。

早期的堤壩潰決研究是以分析潰壩水流運動的理論解為主，1871年Saint-Venant 提出Saint-Venant 方程組后，J.J.Stoker[6]結(jié)合潰壩連續(xù)波和不連續(xù)波研究得到了潰壩波的Stoker 解；林秉南等[7]在假定河槽具有拋物線形斷面且在平底、無阻力的簡化條件下得出了有限長棱柱體水庫的潰壩波對稱解；謝任之[8]概括了現(xiàn)行不同的方法，較系統(tǒng)地提出了對于水庫潰壩流量計算所歸納出的統(tǒng)一公式，給出了相關的表格；進入20世紀80年代后，數(shù)值模擬方面開始逐步廣泛應用，R.Garcia等[9]將數(shù)學模型的數(shù)值解與基于Saint-Venant 方程組的一維無摩擦潰壩解析解進行了對比驗證；R.J.Fennema[10]等采用不連續(xù)的邊界條件模擬二維漸變流；C.Biscarini[11]等提出數(shù)值模擬潰壩水流自由表面流的淺水解決方案；魏文禮等[12]利用 Mac Cormack 預測-校正技術的隱式數(shù)值格式求解二維淺水方程，在此基礎上建立模擬大壩潰決瞬間的洪水演進過程數(shù)學模型，預測矩形河道下游有多個障礙物時大壩局部潰決瞬間的洪水演進過程；冶運濤等[13]基于修正控制方程的復雜邊界潰壩水流數(shù)值模擬，較好模擬復雜邊界處水流特性；程永光[14]歸納總結(jié)了計算流體動力學最新理論知識，通過計算機和數(shù)值方法求解流體力學的控制方程，對流體力學問題進行模擬和分析。

在堤壩潰決后若能以決口現(xiàn)場的水力學特性參數(shù)為依據(jù)迅速有效地制定搶險措施進行封堵，將會大大減少受災區(qū)人民的生命財產(chǎn)損失，本文以某流域河流堤防為例，基于FLOW-3D建立了不同決口口門形狀(U形、V形、T1形和T2形)的堤坊決口水流三維數(shù)值模型，計算并分析了4種口門形狀下決口附近的水位場及流速場分布，所建的數(shù)值模型和模擬方法可為制定最佳堵口方案提供技術參考。

1 基本理論

FLOW-3D為專業(yè)商業(yè)軟件，采用獨創(chuàng)的FAVOR法及真實的Tru-VOF法，具有離散完整的N-S方程，能很好地模擬真實世界中自由液面的流動現(xiàn)象并準確計算其流場。本次數(shù)值模擬選擇非恒定、VOF及RNG湍流模型[15-16]。

1.1 基本方程

基本方程包括連續(xù)性方程和動量方程。

連續(xù)性方程：?(uAx)/?x+?(vAy)/?y+?(wAz)/?z=0

(1)

動量方程：

(2)

(3)

(4)

式中：u，v，w分別為x，y，z方向上的流速分量；Ax，Ay，Az為3個方向上可流動的面積分數(shù)；Gx，Gy，Gz為3個方向上的重力和非慣性力加速度；fx，fy，fz為3個方向上的黏滯力加速度；VF為可流動的體積分數(shù)；ρ為流體密度；p為作用在流體微元上的壓強；t為時間。

1.2 RNG k-ε模型

RNG k-ε模型與標準k-ε模型相比，小尺度運動能從控制方程中去除，可以更好地處理流線彎曲較大及高應變率的流動。

(5)

(6)

1.3 自由表面處理技術

FLOW-3D軟件對自由表面追蹤采用VOF法。其原理[16]是模型中互不相融的兩種或多種流體的任何相位，都由與其相對應的相關變量來描述，計算域整體的流動由所有相位的流動情況相加構(gòu)成。任一單元的變量和特性取決于每一相的面積和體積分量值及其運動情況。

VOF運動學方程：

(7)

流體體積函數(shù)F=F(x,y,z,t)代表計算區(qū)域內(nèi)流體的體積占計算區(qū)域的相對比例。F=1表示該單元完全被流體充滿；F=0表示該單元完全被氣體充滿；0

2 決口水流數(shù)值模型算例模擬驗證

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

為模擬驗證潰壩波的洪水演進過程，選取文獻[9-10]中的局部瞬間潰決水位場計算算例進行驗證。數(shù)值模型尺寸設為200 m×200 m×10 m(圖1)，劃分網(wǎng)格后的單元尺寸為2 m×2 m×1 m，總共1.0×105個單元。模擬的初始狀態(tài)為壩上游水深10 m，下游水深5 m，壩中間缺口處不設置水體，底坡平順且阻尼忽略不計。邊界條件：設重力加速度為-9.8 m/s2，模型頂部設1個標準大氣壓，底部設為無滑移固壁邊界，左右兩側(cè)皆設為固壁邊界。求解設置：求解時間設為10 s，時間步長設為0.02 s。

在模擬時間為7.2 s時，潰壩波波前到達了河道左岸且在河道下游中心發(fā)展完全，因此，國內(nèi)外學者都對其模型在7.2 s 時的水流計算結(jié)果進行分析比較，本文也取7.2 s時的水流計算結(jié)果(水深)與其他文獻相比較。本算例所提取水深計算結(jié)果的截面A-A和B-B及點P1，P2的具體位置見圖1[11]。將模擬計算結(jié)果與文獻[11]中的二維和三維模擬及R.J.Fennema等[10]的散點結(jié)果進行對比(見圖2)。

圖2 水深計算結(jié)果與文獻[11]計算結(jié)果對比Fig.2 Comparison between calculated results of depth given by paper and reference [11]

從圖2可見，計算的水深數(shù)值大小及其變化規(guī)律與文獻[11]中的二維及三維模擬結(jié)果基本吻合。所以，所建數(shù)值模型的計算結(jié)果與其他學者對該算例的計算結(jié)果吻合較好，驗證了所建數(shù)值模型的可靠性以及對決口水流進行三維模擬的可行性。

3 堤防決口三維水流數(shù)值模型

堤防決口一般與天然河床水流方向垂直或斜交，決口呈突發(fā)性，除人工炸堤外，堤防決口一般都不是瞬間潰決，決口都有不同的口門擴寬過程，擴寬至一定程度后會趨于穩(wěn)定?？陂T的擴寬受多種因素影響，擴寬過程十分復雜，為簡化計算，對口門擴寬穩(wěn)定后的堤防決口進行數(shù)值模擬計算，即采用瞬間潰決模式對不同口門形狀下的決口水流進行模擬并分析其結(jié)果。以某實際流域河流堤防進行決口水流模擬建模，然后計算分析不同口門形狀下的堤防決口水流的水力學特性。

圖3 河道及堤防橫斷面(單位：m)Fig.3 Cross-sections of river channel and dike (unit:m)

模型河道內(nèi)的模擬區(qū)域為長700 m，寬100 m(河底)的長直河道，模型高度10 m，河道坡降為0.3‰，河道糙率為0.025，決口高度8 m，距離下游河道200 m，河道左岸堤防高10 m，堤頂寬4 m，兩側(cè)邊坡分別為1∶1.05和1∶0.85，河道及堤防典型斷面如圖3。此外，計算模型還包括決口處50 m內(nèi)的平地段，綜合考慮模型的尺寸和精度，將模型劃分為3個網(wǎng)格塊，第1個網(wǎng)格塊為河道，網(wǎng)格劃分為2 m×2 m×1 m；決口及決口出口50 m為另外2個網(wǎng)格塊，對這2個網(wǎng)格塊進行加密，將網(wǎng)格設為1 m×1 m×1 m。邊界條件：物理邊界條件設重力加速度為-9.8 m/s2，上游邊界設流量為800 m3/s，水位為6 m，下游邊界設為自由出流；模型頂部設1個標準大氣壓，模型底部和模型左右兩側(cè)皆設定為固壁邊界。初始條件：給定河道初始水位為5 m。求解設置：求解時間設為600 s，時間步長設為2 s。

4 不同決口口門形狀水流水力學特性的數(shù)值模擬

4.1 不同決口口門形狀特征

實際決口的口門形狀往往具有任意性，為簡化計算，選擇U形、V型、T1形和T2形這4種典型的決口口門形狀進行數(shù)值模擬。不同口門形狀決口模型的頂部長20 m，寬30 m，高8 m，具體形狀及尺寸如圖4所示。經(jīng)過計算得到不同口門形狀下各決口形狀水流在t=523 s后達到穩(wěn)定，選取t=523 s時刻的計算數(shù)據(jù)進行分析。

圖4 不同形狀決口口門(單位：m )Fig.4 Different shapes of dike breachs (unit：m)

圖5 決口處平面截面示意(單位：m )Fig.5 Dike breaking sections (unit：m)

4.2 水位場分布規(guī)律

選擇決口水流C-C，D-D，E-E，F(xiàn)-F等截面計算結(jié)果進行分析，各截面位置如圖5所示。計算結(jié)果表明，總體上不同口門形狀的決口處水流的水深分布規(guī)律基本相似，水流在進入決口后發(fā)生急速的側(cè)向收縮，在沿決口水流方向(y方向)，水流距離堤腳處越遠，水位越低；在決口寬度方向(x方向)，兩側(cè)堤腳處特別是靠近河道上游來水側(cè)水深較低，而口門中部位置相對較高，水面線在垂直于決口中軸線方向上呈兩邊低、中間高的拱形，水流在決口兩側(cè)堤腳處產(chǎn)生水跌，形成卷吸；各不同決口口門形狀的水深計算結(jié)果見圖6和7。

圖6 不同形狀決口水深計算結(jié)果三維圖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等斷面水面線沿決口寬度方向的分布。不同口門形狀的決口水流水深沿水流方向變化規(guī)律基本一致，均沿水流方向水面迅速降低。在決口進口D-D斷面，4種口門形狀的斷面水深整體較高，水深總體上呈兩邊低中間高的分布特征，且決口右側(cè)水深較左側(cè)水深高；各口門形狀在該斷面水深最大值由大到小分別為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形決口口門寬度最小，故該形狀決口的水面線最窄，最高水位出現(xiàn)在接近決口中軸線附近；其他形狀決口在該兩個斷面的最高水位多數(shù)出現(xiàn)在決口右側(cè)位置，且隨決口斷面面積增大其最高水位越靠近決口右邊界處。

4.3 流速場分布規(guī)律

各口門形狀的決口流速計算結(jié)果見圖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種口門形狀的決口流速分布特征基本相似，在決口處沿水流方向，距堤腳越遠，流速越大；在決口進口橫斷面，最大流速出現(xiàn)在決口左側(cè)，由于部分上游水流發(fā)生轉(zhuǎn)向，決口出口處又無水流的頂托作用，因此流速劇增；在決口出口斷面，水流向四周平地擴散，水面跌落較大，產(chǎn)生較大流速。

決口處斷面流速分布見圖9。

圖9 不同斷面流速分布Fig.9 Velocity distributions along different sections

由圖9可見，C-C斷面上，流速沿水流方向增大；在決口寬度方向，4種口門決口處的流速分布基本相似，均呈不對稱拱形分布。在進口D-D斷面，最大流速均出現(xiàn)在決口中軸線左側(cè)，U形決口在該斷面流速最大，最大流速4.40 m/s，出現(xiàn)位置為x=6.5 m，V形決口在該斷面流速最大值中最小，最大流速為2.45 m/s，出現(xiàn)位置為x=12.5 m。決口E-E和F-F斷面分別靠近下游和出口位置，過流斷面有所減小，流速較進口斷面大，流速沿寬度方向的分布呈中間大兩側(cè)小的拱形分布；對于E-E斷面，除V形決口外，其他口門形狀決口最大流速出現(xiàn)在決口左側(cè)位置；對于決口出口的F-F斷面，流速在3個斷面中最大，各口門形狀流速最大值均在決口中軸線右側(cè)。

5 結(jié) 語

計算分析了在上游流量和水頭及口門寬度一定的條件下不同口門形狀的決口處水位場和流速場分布特征?？傮w上不同口門形狀的決口處水深、流速分布基本相似，數(shù)值上有所差別。各斷面中水深最大值為T1形決口在x=18.5 m處，水深為3.87 m；決口進口斷面最大流速均出現(xiàn)在決口中軸線左側(cè)，最大值為U形決口在x=6.5 m處，流速為4.40 m/s。對決口封堵而言，決口處的流速分布對科學制定決口封堵方案及選擇封堵料粒徑尤為重要，從上述各決口口門形狀的流速分布特征來看，決口口門左側(cè)流速較大，極易造成該側(cè)堤防堤腳的沖刷從而擴大決口寬度，因此在封堵決口時，應及時對該側(cè)決口堤腳進行裹頭處理；在選擇封堵料粒徑時，左側(cè)口門處的封堵粒徑要比右側(cè)口門封堵粒徑大才可保證該側(cè)封堵料的拋投穩(wěn)定。本文所建的決口水流水力學特性模擬的三維數(shù)值模型和模擬方法可為堤防決口的封堵提供技術支持和參考。

[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].湖北水力發(fā)電，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]曲紅玲.河道潰堤與潰堤波的一、二維耦合計算數(shù)值模擬[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]魏文禮,沈永明,孫廣才,等.二維潰壩洪水波演進的數(shù)值模擬[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]冶運濤,梁犁麗,張光輝,等.基于修正控制方程的復雜邊界潰壩水流數(shù)值模擬[J].水力發(fā)電學報,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]張婷.波浪的三維數(shù)值模擬及其應用[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

李火坤，鄧冰梅，曾秀娟，等.堤防決口水流水力學特性的三維數(shù)值模擬[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) ；廣東省水利科技創(chuàng)新基金資助項目(2014-08)

李火坤(1981—)，男，湖南長沙人，博士，副教授，主要從事水工水力學及泄流結(jié)構(gòu)動力檢測方面的研究工作。E-mail:lihuokun@126.com

TV871； TV122+.4

A

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