王 兵，張鳳德，賀昌海

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

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

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

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

1 研究河道概況

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

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

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

2 數(shù)值計(jì)算

2.1 控制方程

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

連續(xù)方程：

(1)

動(dòng)量方程：

(2)

紊動(dòng)能k方程：

(3)

紊動(dòng)能耗散率ε方程：

(4)

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

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

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

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

圖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è)為壓力邊界，用下游水位控制，下游水位通過(guò)試算確定；左岸Xmin邊界與流體不接觸，不影響流態(tài)，采用缺省設(shè)定為對(duì)稱邊界；潰口堤后Xmax設(shè)為壓力邊界，根據(jù)與百年一遇洪水位(潰口處河道斷面)的水位落差ΔZ1控制；模型底部Zmin設(shè)為固壁邊界；模型頂部Zmax設(shè)為壓力邊界，P=0(一個(gè)大氣壓)。Mesh②除頂部與底部外，均取對(duì)稱邊界。

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

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

2.4 計(jì)算工況

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

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

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

3.1 不同口門(mén)寬度下流速場(chǎng)

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

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

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

圖4 口門(mén)區(qū)域流速場(chǎng)Fig.4 Velocity field of the entrance area

3.2 不同口門(mén)寬度下水位場(chǎng)

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

圖5 口門(mén)中軸線水面線Fig.5 Water surface profile of the central axis of the entrance

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

圖6 口門(mén)區(qū)域水位場(chǎng)Fig.6 Free surface elevation of the entrance area

3.3 局部大流速分布特征

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

圖7 口門(mén)斷面流速分布圖Fig.7 Velocity distribution diagram of breach section

3.4 潰口分流量和口門(mén)斷面平均流速

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

4 結(jié) 語(yǔ)

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

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

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

□

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

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

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

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

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

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

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

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

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

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

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

[12] 劉成林，陳宇豪．基于Flow-3D的水平射流沖刷泥沙數(shù)值模擬[J]．人民長(zhǎng)江，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é)報(bào)(自然科學(xué)與工程技術(shù)版)，2016，49(4)：422-428．