基于Flow-3D的堤防堵口水力條件數(shù)值模擬

王 兵，張鳳德，賀昌海

(1. 武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 黑龍江省三江工程建設(shè)管理局，哈爾濱 150081)

堤壩潰決水流一般兼有急流和緩流、涌波和稀疏波等流態(tài)，決口流量過程常是峰高量大，形狀陡峭，急緩流轉(zhuǎn)換，水面梯度較大，在口門區(qū)域常存在間斷[1,2]。潰堤問題與潰壩問題相比最大的區(qū)別在于河道的側(cè)向水流運動，與側(cè)堰水流運動非常相似[3]。目前針對河道堤防潰決水流的相關(guān)數(shù)值模擬研究，主要為二維平面數(shù)值模擬[4-6]。

為了盡量減少堤防潰決后造成的損失，結(jié)合實際地形地質(zhì)條件[7-9]和相應(yīng)口門區(qū)域水力條件及變化規(guī)律，制定科學(xué)合理的堵口搶險預(yù)案在當(dāng)前時代特征下顯得尤為重要。在傳統(tǒng)物理模型試驗存在周期長、投資大等缺點的背景下，近年來Flow-3D等CFD軟件得到了越來越廣泛的運用，包括溢流壩泄流[10]、消能池水躍[11]、泥沙輸運[12]、長廊式氣墊調(diào)壓室的水-氣兩相流[13]等三維數(shù)值模擬計算。眾多前人的研究成果表明，三維水流數(shù)值模擬計算結(jié)果與物理模型試驗結(jié)果是相對一致的[14]，可用于指導(dǎo)和解決實際工程問題。

本文利用Flow-3D軟件，針對某河道砂性堤防潰決水流運動進(jìn)行三維數(shù)值模擬，研究超標(biāo)準(zhǔn)洪水條件下口門區(qū)域水力條件，其成果可為研究極端洪水條件下該河道砂性堤防潰口快速封堵技術(shù)提供科學(xué)依據(jù)。

1 研究河道概況

某河道兩岸地形地貌為低漫灘階地，灘地高程多在41～44 m之間，而島面高程多為41～42 m。河床與河岸的土質(zhì)以細(xì)砂及壤土為主，抗沖刷能力差。平槽水位時河道水深達(dá)10 m以上，水面寬闊，水深流急。該河道堤防防洪標(biāo)準(zhǔn)為50 a一遇，堤防形式以砂基砂堤或砂基土堤為主，部分堤段受地形地質(zhì)等因素的影響為險工弱段。

研究選取某個典型堤防險段及其對應(yīng)的河道實際地形，對百年一遇洪水條件下的堤防潰口水流運動進(jìn)行數(shù)值模擬。該典型險段設(shè)計斷面如圖1，堤頂超高2.2 m，堤頂寬8 m，兩側(cè)邊坡均為1∶4.0。

圖1 典型險段設(shè)計斷面(單位：cm)Fig.1 Design section of typical risk dike

2 數(shù)值計算

2.1 控制方程

本文研究的具有實際復(fù)雜地形的堤防潰口水流問題，因地勢起伏多變及河道的側(cè)向水流運動特性，潰口區(qū)域水流流態(tài)復(fù)雜、水面梯度大，因此采用RNGk-ε紊流模型?？刂品匠倘缦拢?/p>

連續(xù)方程：

(1)

動量方程：

(2)

紊動能k方程：

(3)

紊動能耗散率ε方程：

(4)

式中：ui是流速張量分量；Ai、gi、fi分別代表三維直角坐標(biāo)方向上可流動的面積分?jǐn)?shù)、重力加速度和黏滯力；VF是可流動的體積分?jǐn)?shù)；ρ是流體密度；p是作用在流體微元上的壓力；k、ε分別是紊動能和紊動能耗散率；σk、σε分別是紊動能和耗散率對應(yīng)的Prandtl數(shù)；μ、μt分別是水的動力黏滯系數(shù)和紊動黏滯系數(shù)；Gk是紊動能k的產(chǎn)生項；C*ε1、Cε2是經(jīng)驗常數(shù)。

采用交錯矩形網(wǎng)格的FAVOR法離散控制方程，采用TruVOF方法處理流體自由表面，采用極小殘差(GMRES)算法求解壓力。

2.2 模型及網(wǎng)格劃分

在原有地形資料的基礎(chǔ)上，利用CATIA軟件進(jìn)行三維河道地形建模，并創(chuàng)建堤防建筑物。將生成的模型實體導(dǎo)入Flow-3D中進(jìn)行數(shù)值計算。為了保證計算的效率、穩(wěn)定性和準(zhǔn)確性，采用均勻的立方體網(wǎng)格，縱橫比為1；并考慮深度方向的計算精度要求相對較高，設(shè)置網(wǎng)格尺寸X∶Y∶Z=2∶2∶1。為更精確的模擬潰口區(qū)域的水力條件，在潰口區(qū)域采用嵌套網(wǎng)格。網(wǎng)格塊①覆蓋河道范圍對應(yīng)堤線樁號從1+000 m至3+150 m處(原型)，立方體網(wǎng)格尺寸為4 m×4 m×2 m；嵌套網(wǎng)格塊②覆蓋范圍為潰口口門堤線前后200 m、上下堤頭各退后50 m所圍成的區(qū)域，潰口中心對應(yīng)堤線樁號2+100 m處，立方體網(wǎng)格尺寸為2 m×2 m×1 m。因此網(wǎng)格總數(shù)隨口門寬度變化而稍有變化，例如口門寬60 m時，網(wǎng)格總數(shù)約379萬，有效網(wǎng)格數(shù)約為229萬。模型及網(wǎng)格劃分如圖2所示，(a)中網(wǎng)格單元擴大六倍以方便顯示。

圖2 模型及網(wǎng)格劃分示意圖Fig.2 Schematic diagram of model and grid division

2.3 邊界條件及初始條件

邊界條件(Mesh①)：上游Ymin設(shè)為流量邊界，根據(jù)所建模型地形斷面數(shù)據(jù)和百年一遇洪水位46.17 m換算而得；下游Ymax設(shè)為壓力邊界，用下游水位控制，下游水位通過試算確定；左岸Xmin邊界與流體不接觸，不影響流態(tài)，采用缺省設(shè)定為對稱邊界；潰口堤后Xmax設(shè)為壓力邊界，根據(jù)與百年一遇洪水位(潰口處河道斷面)的水位落差ΔZ1控制；模型底部Zmin設(shè)為固壁邊界；模型頂部Zmax設(shè)為壓力邊界，P=0(一個大氣壓)。Mesh②除頂部與底部外，均取對稱邊界。

初始條件：以下游控制水位作為計算域內(nèi)的初始水位。

流量監(jiān)測：Baffles 在FLOW-3D中是沒有厚度的孔隙孔板，用以控制或引導(dǎo)水流。實際工作中常用于測量通過某斷面的水流流量與計算通過的顆粒數(shù)量，且完全不影響水流運動。因此設(shè)置Baffle1、Baffle2和Baffle3分別監(jiān)測模型上游進(jìn)口流量Q1、下游出口流量Q2和潰口處分流量QK，如圖2所示。

2.4 計算工況

研究主要利用FLOW-3D分別對不同口門寬度及堤防內(nèi)側(cè)不同水位條件作用下的潰口水流運動進(jìn)行數(shù)值模擬，得到不同工況條件下口門區(qū)域的水位、水位場分布、流速及流速場分布等水力學(xué)特性，進(jìn)而分析其對口門區(qū)域水力條件的影響程度。計算工況如表1。

表1 數(shù)值模擬計算工況Tab.1 Calculation conditions of numerical simulation

3 數(shù)值計算結(jié)果分析

3.1 不同口門寬度下流速場

圖3為沿堤防軸線口門斷面的沿水深平均流速分布圖，在不同口門寬度下其流速分布形態(tài)基本一致?？陂T斷面上，靠近兩堤頭流速急劇減小，中間段流速分布較堤外斷面更為均勻，但下堤頭前出現(xiàn)局部流速峰值，而數(shù)值隨口門束窄先增大后減小。圖3(b)與圖3(a)相比可知，隨著ΔZ1減小，也即堤后水位上漲，口門寬度相同條件下斷面流速明顯減小。

圖3 口門斷面沿水深平均流速Fig.3 Depth-averaged velocity of the entrance section

堤外流速峰值分布于口門偏上游區(qū)域，隨著口門的束窄流速變化坡度增大。堤內(nèi)流速集中于口門中間區(qū)域，橫向擴散現(xiàn)象明顯。例如口門寬60 m時口門區(qū)域流速場，如圖4。

圖4 口門區(qū)域流速場Fig.4 Velocity field of the entrance area

3.2 不同口門寬度下水位場

圖5為沿口門中軸線斷面水面線分布圖，潰口水流水面梯度隨口門寬度的減小而增大。堤線前水位隨口門寬度減小而升高，水流經(jīng)過口門過程中出現(xiàn)垂向收縮的現(xiàn)象。

圖5 口門中軸線水面線Fig.5 Water surface profile of the central axis of the entrance

在上、下堤頭前，堤角附近的卷吸現(xiàn)象明顯，對口門區(qū)域橫向水位變幅影響較大。例如口門寬60 m時口門區(qū)域水位場，如圖6。

圖6 口門區(qū)域水位場Fig.6 Free surface elevation of the entrance area

3.3 局部大流速分布特征

根據(jù)數(shù)模計算結(jié)果，堵口進(jìn)占過程中口門區(qū)域存在局部大流速分布。60 m龍口寬度時口門斷面流速分布如圖7，下堤頭附近流速明顯高于上堤頭，工況7條件下最大流速達(dá)8.27 m/s，工況8條件下最大流速達(dá)6.74 m/s，位置向下堤頭前靠攏；同時口門底部流速較大。因此合龍前必須分別針對下堤頭和口門底部的沖刷情況采取相應(yīng)的裹護(hù)措施。

圖7 口門斷面流速分布圖Fig.7 Velocity distribution diagram of breach section

3.4 潰口分流量和口門斷面平均流速

根據(jù)Baffles監(jiān)測結(jié)果，堤防潰口口門寬度與堤后邊界條件對潰口分流量的影響十分顯著。各工況條件下潰口分流量與口門斷面平均流速如表2。

4 結(jié) 語

研究選取某河道堤防典型險段實際地形和現(xiàn)有防洪設(shè)計標(biāo)準(zhǔn)(P=2%)的堤防建筑物，采用CATIA軟件作為三維建模平臺，并基于FLOW-3D軟件對百年一遇洪水條件下、不同口門寬度及堤防內(nèi)側(cè)不同水位條件作用下的潰口水流進(jìn)行數(shù)值模擬，從而得到14種組合工況條件下口門區(qū)域的水位、水位場分布、流速及流速場分布等堤防堵口水力學(xué)特性。計算結(jié)果表明，潰口內(nèi)外水位差和口門寬度均對潰口分流量、口門斷面平均流速有顯著的影響，應(yīng)當(dāng)據(jù)此選擇合適的堵口時機；隨著口門束窄，潰口前水位逐漸抬高，口門斷面流速先增后減，堵口進(jìn)占過程中因根據(jù)流速變換拋投材料及粒徑；合龍過程中應(yīng)適時監(jiān)測口門底部和下堤頭前的沖刷現(xiàn)象，并采取相應(yīng)的防護(hù)措施；為保證堵口截流順利實施，應(yīng)當(dāng)具備足夠的拋投強度等。

表2 潰口分流量和口門斷面平均流速Tab.2 Flow rate and average velocity of breach section

注：堤身內(nèi)外實際落差ΔZ2以口門中軸線上距離原堤軸線100 m的兩點處水位差計。

□

[1] 譚維炎．計算淺水動力學(xué)[M]．北京：清華大學(xué)出版社，1998．

[2] 張大偉．堤壩潰決水流數(shù)學(xué)模型及其應(yīng)用研究[D]．北京：清華大學(xué)，2008.

[3] 張修忠．堤防潰決過程的數(shù)值模擬[D]．北京：清華大學(xué)，2003.

[4] 苑希民,田福昌,王麗娜．漫潰堤洪水聯(lián)算全二維水動力模型及應(yīng)用[J].水科學(xué)進(jìn)展，2015，26(1)：83-90.

[5] 賀治國.考慮水沙相互作用的河堤潰決二維數(shù)值模型[J].水動力學(xué)研究與進(jìn)展A輯，2010，25(2)：147-154.

[6] 陳 珺，張小峰，談廣鳴，等.考慮潰口展寬的潰堤水流泥沙數(shù)值模擬[J].水動力學(xué)研究與進(jìn)展A輯，2007，22(5)：647-653.

[7] 楊佩國，楊勤業(yè)，吳紹洪,等．基于數(shù)值模擬的黃河下游不同情景潰堤洪水特性[J]．地理研究，2007，26(2)：328-336．

[8] 張大偉，張 超，王興奎．具有實際地形的潰堤水流數(shù)值模擬[J]．清華大學(xué)學(xué)報(自然科學(xué)版)，2007，47(12)：2 127-2 130．

[9] 徐國賓，孟慶林，苑希民．含沙潰堤洪水?dāng)?shù)值模擬分析[J]．天津大學(xué)學(xué)報(自然科學(xué)與工程技術(shù)版)，2016，49(10)：1 008-1 015．

[10] 戚 藍(lán)，陳 輝，費文才，等．基于精確河道地形的溢流壩泄流三維數(shù)值模擬[J]．水力發(fā)電學(xué)報，2016，35(10)：12-20．

[11] 王月華，包中進(jìn)，王 斌．基于Flow-3D軟件的消能池三維水流數(shù)值模擬[J]．武漢大學(xué)學(xué)報(工學(xué)版)，2012，45(4)：454-457．

[12] 劉成林，陳宇豪．基于Flow-3D的水平射流沖刷泥沙數(shù)值模擬[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建模的導(dǎo)流工程三維數(shù)值模擬研究[J]．天津大學(xué)學(xué)報(自然科學(xué)與工程技術(shù)版)，2016，49(4)：422-428．