曾秀娟，謝 濤

(1.江西省上饒市水利電力勘測設計院，江西 上饒 334000；2.江西省贛州市水利電力勘測設計研究院，江西 贛州 341000)

堤防是世界上最早廣為采用的一種重要防洪工程，是防洪體系中重要的擋水建筑物。堤防給人們帶來了巨大的經濟效益和社會效益，但因漫溢、沖刷、滲透、凌汛、超標準洪水等原因造成堤防決口時，也會帶來巨大的損失。在眾多洪水險情中，堤防決口是損失最嚴重、影響面最大、搶險最艱難的[1]。在堤防決口后，若能迅速制定有效的相關封堵決口措施，將大大減少決口造成的損失。通過對堤防決口進行不同工況下的模擬、研究，得到堤防決口處的流速、水面線等水力學特性的分布和變化規(guī)律，從而針對性的制定有效的封堵方法，具有重要的意義。近年來，已有部分堤防決口的數值模擬研究成果：李火坤，鄧冰梅等[2]對不同決口口門形狀下的堤防決口水流建立了三維數值模型；李火坤、曾智超等[3]對堤防決口采用立堵法、平堵法的封堵方法進行數值模擬，得到封堵過程中決口附近的水位場和流速場分布規(guī)律；羅娜、劉成林等[4]對不同型式下堤防決口的沖刷過程進行了數值模擬。本文采用FLOW-3D軟件以某實際河道堤防為模型，建立了不同水頭(6 m、7 m、8 m、9 m)、不同流量(400 m3/s、600 m3/s、800 m3/s、1 000 m3/s和1 500 m3/s)作用下的堤防決口三維數值模擬，得到4種不同水頭、5種不同流量下的決口處水流的水力學特性，為擬定有效的決口封堵方案提供科學依據和技術參考。

1 FLOW-3D軟件基本介紹

1.1 FLOW-3D軟件

FLOW-3D是一款計算流體力學的軟件，近年來在水利方面已有較廣泛的應用[2-6]，它采用獨創(chuàng)的FAVOR(結構化矩形網格)方法及真實的Tru-VOF(可視化)方法[7]，可以模擬真實世界中物理模型的流動現象及準確計算出各種流場性質，特別是自由液面的流動[2]。本文計算模型選擇：非恒定模型，VOF模型和RNG湍流模型。

1.2 基本方程

基本方程主要包括連續(xù)性方程和動量方程[2-4]。

連續(xù)方程式：

(1)

動量方程式：

(2)

(3)

(4)

式中：RDIF為密度擴散；U=(u，v，w)為液體速度，m/s；P為壓力，kPa；G為重力和非慣性力加速度，m/s2；τ為粘性應力張量；K為阻力、拖曳力，kN；RSOR為在流速為零的時刻由于大量的注射劑引起的加速度，m/s2；F為其它力(包括表面張力、電場力、動力來源、粒子及用戶自定義的力)，kN。

1.3 自由表面處理技術

FLOW-3D軟件中自由液面是指具有液氣共存的介面模型，在氣體中的壓力梯度及剪應力是忽略的，通常應用于密度比高達1 000∶1的兩種流體。FLOW-3D中自由液面的簡化模型：(1)忽略氣體流動，氣體僅有一個正向壓力作用于液體表面上；(2)自由表面上的氣體以壓力邊界取代。FLOW-3D軟件采用VOF 法對自由表面進行追蹤。 對于水氣二相流而言，其基本思想是：定義函數αω(x，y，z，t)和αa(x，y，z，t)分別代表計算區(qū)域內水和氣占計算區(qū)域的體積的相對比例。αω=1，表示該單元完全被水充滿；αω=0，表示該單元完全被氣充滿；αω<1，表示該單元部分是水，部分是氣，有自由面。

2 模擬驗證

本文選用有較大模擬結果的文獻[8]中的一個有矩形障礙物的潰壩算例進行驗證。模型尺寸為500 m×300 m×10 m，網格劃分為2 m×2 m×1 m，網格總數為3.75×105個。算例給定初始條件為上下游初始水位分別為10 m、5 m，壩中間缺口處水位為0。物理邊界條件：重力Z方向給定一個重力加速度，x、y方向不做設置；模型的頂部和上下游邊界都設為壓力邊界條件，左右兩側和底部邊界都設定為固壁邊界。求解時間設為40 s，時間間隔設為0.5 s。模型平面尺寸見圖1所示。

圖2從左往右依次是FLUENT[8]、Wang[9]、FLOW-3D的計算時間分別為12 s、18 s、32 s的速度場計算結果。由圖可知，FLOW-3D的模擬計算成果與其他兩個文獻的結果非常相似：在12 s時刻，洪水都未到達矩形障礙物；在18 s時刻，洪水剛到達障礙物；在32 s時刻，洪水波到達下游邊界并對矩形障礙物形成繞流，在模型的下游范圍內水面相互交叉重疊，從圖中可以清楚的看到波的反射和波與波之間的相互作用。通過以上分析可知，用FLOW-3D所建的數值模型計算成果與該算例的其他成果比較吻合，驗證了Flow-3D軟件用于堤防決口三維數值模擬的可行性。

圖1 算例模型平面尺寸

3 模型建立

本文采用某實際河道堤防簡化建模，模型尺寸為700 m×160 m×10 m(長×寬×高)，模型包括尺寸為700 m×100 m(河底)×10 m(長×寬×高)的河道、河道兩岸堤防與堤腳處等高的尺寸為700 m×50 m×10 m(長×寬×高)的平地；模擬河道坡降和糙率分別設為0.3‰和0.025，河道堤防高10 m，堤頂寬4 m，兩側邊坡分別為1∶1.05、1∶0.85，設置決口處距離下游河道200 m[9]；河道堤防斷面及Flow-3D模型三維圖見圖3。

河道模型網格劃分為2 m×2 m×1 m，在決口50 m×50 m平面范圍內進行網格局部加密，加密網格為1 m×1 m×1 m。模型初始條件為設置河道水位為5 m；物理邊界條件：重力Z方向給定一個重力加速度，x、y方向不做設置；不同水頭工況下模擬設置初始流量為800 m3/s，水位隨工況變，不同流量工況下模擬設置初始水頭6 m，流量隨工況變；下游河道設為自由出流；模型頂部設置壓力為一個標準大氣壓壓力邊界條件，模型底部和左右兩側皆設定為固壁邊界。求解時間設為600 s，時間步長設為2 s。決口口門三維結構尺寸及平面尺寸見圖4，圖中x方向為河道水流方向(模型中坐標為470～500)，y方向為決口水流方向(模型中坐標為100～120)，選擇圖4中的A-A和B-B截面的水深和流速計算結果進行分析。

圖2 FLUENT、Wang、FLOW-3D計算的速度場結果(時間分別為12 s、18 s、32 s)

圖3 河道堤防斷面及Flow-3D模型三維圖

圖4 模型決口口門平面尺寸及三維結構圖

4 不同水頭作用下堤防決口的水力學特性分析

本文分別對河道上游水頭為6 m、7 m、8 m、9 m的條件下T形決口(圖4)為例進行堤防決口水力學特性的數值模擬。表1為不同水頭作用下決口處水流流態(tài)達到穩(wěn)定所需的時間，由表可知，在給定河道初始水位為5 m、上游來水流量不變的條件下，決口水流模擬達到穩(wěn)定的時間不會隨上游水頭的增大而發(fā)生改變。本文主要對不同水頭作用下水流流態(tài)達到穩(wěn)定后的決口處平面流速場分布規(guī)律及A、B各斷面的水面線和流速水力數據進行分析。

表1 不同水頭作用下數值模擬穩(wěn)定時間

4.1 流速場分布

圖5是不同水頭作用下決口平面三維流速場。從圖中可以看出，不同水頭作用下的決口流速場分布幾乎完全一致，水流在決口進口左右兩側產生側向收縮，水流產生明顯的縱向跌落；因模型設置決口出口為無任何阻擋物的平地，水流沿決口口門出口在重力勢能作用下向四周擴散，故沿決口水流方向，流速逐漸增大；在垂直于水流方向上，因口門形狀比較規(guī)整，在堤腳處產生了水跌，水流紊動劇烈，流速較大。從圖中可以明顯看到決口平面處的較大流速主要分布在決口進口上游左側發(fā)生水跌處(圖中紅色部分)，決口進口右側河道內水流流速幾乎為0(圖中藍色部分)。在進行決口封堵時，當決口來水側水流流速大，水流運動紊亂，應考慮在該側進行裹頭處理，從決口下游側采用單邊進占封堵，直至合龍。

圖5 不同水頭下T形決口三維流速場

4.2 A截面水面線和流速結果分析

圖6和圖7分別是不同水頭作用下A-A斷面的水面線和流速。從圖中可以看到，不同水頭作用下的水面線沿水流方向產生縱向跌落，水位沿水流方向距離堤腳處越遠，水位越低，水位逐漸降低，水面線為正比降，四種水頭下的水面線基本平行；水流沿水流方向距離堤腳處越遠，流速越大。在106.5 m≤y≤116.5 m之間，各不同水頭作用下流速分布是6 m水頭<7 m水頭<8 m水頭<9 m水頭；在y≤106.5 m、116.5 m≤y范圍內，各不同水頭作用下流速分布規(guī)律不明顯，流速線有部分交叉，較大流速分布在決口下游出口邊界處。隨著作用水頭的增加，決口中軸線A-A斷面處的水位和流速隨之加大。

圖6 不同水頭下A-A斷面水面線

圖7 不同水頭下A-A斷面流速

4.3 B截面水面線和流速結果分析

圖8和圖9分別是B-B斷面的水面線和流速。B-B斷面的水面線呈現左邊峰值較低、右邊峰值稍高的高拱形不對稱分布，和流速分布正好相反，各不同水頭作用下潰口的水面線、流速變化規(guī)律一致，且數值都比較接近，水面線互有一定交叉；從表2中和圖8、9中可以看出，8 m水頭作用下的B截面斷面流速均稍大于9 m水頭作用下的B截面斷面流速，這是因為在上游流量一定的情況下，進入決口的流量相差不大，而T形決口隨著決口處水面的增高其過水斷面面積增大，在流量相差不大的情況下，過水面積增大，流速減小。從B-B斷面最大水位和最大流速分布(表2：流速單位m/s，水深單位m)可以看出，決口較大的流速分布在決口進口上游側，隨著作用水頭的增加，堤防中軸線B-B斷面的水位和流速隨之加大。

表2 不同水頭下B-B斷面最大水位和最大流速分布

圖8 不同水頭下B-B斷面水面線

圖9 不同水頭下B-B斷面流速

5 不同流量作用下堤防決口的水力學特性分析

河道上游不同流量條件對決口處的水力學特性參數也會產生不同的影響，本文分別對河道上游流量為400 m3/s、600 m3/s、800 m3/s、1 000 m3/s和1 500 m3/s的條件下進行堤防決口水力學特性的數值模擬。表3為不同流量下決口處水流流態(tài)達到穩(wěn)定所需的時間。由表3可知，在給定河道初始水位為5 m、上游水頭不變的條件下，決口水流模擬達到穩(wěn)定的時間隨上游來流量的增大而減小。因不同流量作用下水流流態(tài)達到穩(wěn)定后的決口處平面流速場分布規(guī)律與不同水頭作用下流速場分布規(guī)律比較類似，這里不再贅述，本節(jié)僅對A、B各斷面(見圖4)的水面線和流速水力數據進行分析。

表3 不同流量作用下數值模擬穩(wěn)定時間

5.1 A截面水面線和流速結果分析

圖10和圖11是不同流量作用下A-A斷面的水面線和平均流速。從該圖可以看到，不同流量作用下的水面線沿水流方向產生縱向跌落，決口進口斷面和出口斷面水位相差較大，各不同流量作用下的水面線基本為線性線段，基本互相平行，決口處水深隨上游來流量的增加而增大。各不同流量作用下決口處各斷面的流速都是沿水流方向逐漸增大，斷面上各點流速隨距堤腳線距離越遠基本成線性增長，流速變化規(guī)律相似，較大流速分布在決口下游出口邊界處。

圖10 不同流量下A-A斷面水面線

圖11 不同流量下A-A斷面流速

5.2 B截面水面線和流速結果分析

圖12和圖13是堤防中軸線斷面B-B斷面的水面線和平均流速。該斷面的水面線呈中間高、兩邊低的高拱形，決口左側水位都比較低，中軸線附近及偏口門右側水位較高；流速呈左邊峰值較高、右邊峰值稍低的城門洞形不對稱分布；各不同流量作用下潰口的水面線和平均流速變化規(guī)律一致。從B-B斷面最大水位和最大流速分布(表4流速單位m/s，水深單位m)可以看出，決口較大的流速分布在決口進口上游側，隨著上游流量的增加，堤防中軸線斷面處的水位和流速隨之加大。在縱向上決口水位隨距堤腳線距離越遠而減小，流速在橫向上從左到右先增大后減小，決口左側邊界附近流速較大，中間部分流速較小。

表4 不同流量下B-B斷面最大水位和最大流速分布

圖12 不同流量下B-B斷面水面線

圖13 不同流量下B-B斷面流速

6 結 論

本文采用Flow-3D對某實際河道的堤防進行決口三維模擬，主要分析了在口門形狀及大小一定的條件下不同水頭、不同流量作用下的決口處水深和流速的分布特征，得到如下結論：(1)水位：各不同水頭、不同流量作用下的水位在橫向上是潰口左右兩側邊界處低、中間高，水流潰口水面線在橫向上呈兩邊低、中間高的拱形。(2)流速：不同水頭、不同流量作用下，潰口流速分布在橫向上呈兩邊低、中間高的拱形；決口較大的流速主要分布在決口進口上游左側發(fā)生水跌處。因此，在進行決口封堵時，應注意預先做好決口上游側的裹頭處理，防止在封堵的過程中因流速過大導致潰口口門擴寬；在選取封堵材料時，也可根據潰口處流速分布來選擇相應的封堵料粒徑。本文基于Flow-3D建立的堤防決口三維模型和模擬方法均可為制定有效的決口封堵方案提供科學依據和技術參考。

[1] 孫蘆忠，趙建均，嚴建國，等.堤防決口的水力學試驗研究[J].人民長江，2003，34(11)：41～42.

[2] 李火坤，鄧冰梅，曾秀娟，等.堤防決口水流水力學特性的三維數值模擬[J].水利水運工程學報，2016(6)：29～36.

[3] 李火坤，曾智超，鄧冰梅，等.堤防決口封堵的水力學特性研究[J].水利水運工程學報，2017(3)：8～15.

[4] 羅娜，劉成林，陳宇豪.堤防決口下游泥沙沖刷的三維數值模擬分析[J].南昌大學學報(工科版)，2016，38(4)：360～365.

[5] 王月華，包中進，王斌.基于FLOW-3D軟件的消能池三維水流數值模擬[J].武漢大學學報(工科版)，2012，45(4)：454～457.

[6] 包中進，王月華.基于FLOW-3D軟件的溢洪道三維水流數值模擬[J].浙江水利科技，2012(2)：5～9.

[7] 曾秀娟.堤防決口封堵水力學數值模型研究[D].南昌大學，2015.

[8] 林長強.基于FLUENT的土石壩逐漸潰壩水流模擬[D].華中科技大學，2011.

[9] J.S.Wang，H.G.Ni，Y.S.He.Finite-Difference TVD Scheme for Computation of Dam-Break Problems.Journal of Hydraulic Engineering，2000，126(4)：253～262.