MIKE軟件技術(shù)在水庫(kù)潰壩洪水?dāng)?shù)值模擬中的應(yīng)用研究

管 靜

(河北省承德水文勘測(cè)研究中心,河北 承德 067000)

0 引 言

水庫(kù)大壩發(fā)生潰壩的概率雖然較低,但由于我國(guó)氣候和地勢(shì)因素,超負(fù)荷的洪水與強(qiáng)震同樣也會(huì)導(dǎo)致潰壩的形成[1-2]。潰壩洪水的流動(dòng)狀況十分復(fù)雜,很多因素都會(huì)影響到潰壩的洪水。目前,對(duì)其模擬主要包含物理模擬、數(shù)學(xué)模擬以及二者結(jié)合的模擬方法[3-4]。其中,數(shù)學(xué)模擬因?yàn)槠潇`活強(qiáng)、資金消耗較少等優(yōu)勢(shì)在國(guó)內(nèi)外得到廣泛的應(yīng)用[5]。高普陽(yáng)等[6]針對(duì)河流下游存在相關(guān)障礙物的潰壩流動(dòng)實(shí)際問(wèn)題,在兩相流動(dòng)模型與有限元算法的基礎(chǔ)上,提出了牛頓流體潰壩流動(dòng)相關(guān)的數(shù)值模擬方法。Maghsoodi R[7]針對(duì)潰壩障礙物對(duì)河流影響的相關(guān)問(wèn)題,在流體體積和標(biāo)準(zhǔn)k-1模擬水面的相關(guān)方法基礎(chǔ)上,對(duì)其進(jìn)行了數(shù)值模擬。呂松峰等[8]針對(duì)三維地形尾礦庫(kù)潰壩的相關(guān)問(wèn)題,在顆粒流離散元法與短矩陣離散元軟件的基礎(chǔ)上,構(gòu)建了實(shí)際云南某尾礦庫(kù)完整的三維地形,并以此對(duì)其下泄演進(jìn)進(jìn)行了數(shù)值模擬。

在此背景下,本文引入土壩失事洪水演算模型(BREACH)與二維數(shù)值工程軟件(MIKE21),并將二者組合構(gòu)建成耦合模型,同時(shí)將其運(yùn)用在河北邢臺(tái)水庫(kù)潰壩的洪水?dāng)?shù)值模擬應(yīng)用中。其目的是為邢臺(tái)水庫(kù)在形成潰壩之前制定合理的防洪避險(xiǎn)方案,并通過(guò)預(yù)測(cè)潰壩對(duì)下游的影響,為下游制定進(jìn)一步的防洪計(jì)劃和轉(zhuǎn)移方案提供理論依據(jù)。此外,通過(guò)將耦合模型應(yīng)用在潰壩洪水?dāng)?shù)值模擬過(guò)程中,制定7種不同的方案,既能全方位分析不同情況的影響,又具備一定的創(chuàng)新性。

1 基于MIKE技術(shù)的水庫(kù)潰壩洪水?dāng)?shù)值模擬分析

1.1 水庫(kù)潰口流量計(jì)算模型研究

在實(shí)際的水庫(kù)潰壩中,不同的大壩類(lèi)型會(huì)導(dǎo)致不同的崩塌形式。根據(jù)壩體的潰決歷史,可以將其分為兩類(lèi),即暫時(shí)性潰決和漸進(jìn)潰決[9-10]。通常來(lái)講,混凝土壩體可以按瞬間潰決形式進(jìn)行計(jì)算,而土石壩往往采用逐步潰決模型進(jìn)行模擬分析。因此,在確定潰壩方式之后,需要選取合適的數(shù)學(xué)模型來(lái)動(dòng)態(tài)分析潰壩的大小和流速,并以此來(lái)模擬潰決下游的洪水演變。目前,用于水庫(kù)潰壩洪水?dāng)?shù)值模擬的數(shù)學(xué)模型有很多,研究依據(jù)模型適應(yīng)能力水平,選擇BREACH模型。

BREACH模型可以對(duì)潰口特性的演變過(guò)程以及潰決后的泄洪過(guò)程進(jìn)行預(yù)測(cè),同時(shí)也可以用于模擬由于漫頂或管涌而造成的潰壩,壩體可以是均勻的,也可以是由兩種材料組成的心壁和外壁[11-12]。BREACH模型的基本原理是將水力學(xué)、泥沙運(yùn)動(dòng)、土力學(xué)原理等因素進(jìn)行綜合考慮,首先假定潰口是長(zhǎng)方形的,然后在潰口處的崩塌和變形中形成一個(gè)梯形,當(dāng)潰口超過(guò)一定深度時(shí),洪水就會(huì)停止沖刷潰口處的土體。潰口不再發(fā)展,而是趨向一個(gè)穩(wěn)定的形態(tài),由此可以預(yù)測(cè)堤壩的潰口大小和潰口的流速。

BREACH模型可以細(xì)分為漫頂潰決和管涌潰決。對(duì)于漫頂潰決來(lái)講,當(dāng)坡面沒(méi)有植被覆蓋時(shí),由于水流的沖刷,在下游坡面上逐漸形成一個(gè)矩形的溝道。此時(shí),河流水流流量計(jì)算公式如下:

Pa=3A0(h-hb)

(1)

式中:Pa為潰口河道的實(shí)際流量;A0為初始矩形河道的瞬間寬度;h為壩前的實(shí)際水位高度;hb為河道底部的實(shí)際高度。

因此,水壩潰口的實(shí)際發(fā)展過(guò)程見(jiàn)圖1。

圖1 水壩潰口的發(fā)展示意圖

由圖1可知,在植被覆蓋的情況下,如果坡面上的水流速度超過(guò)允許速度,下游的坡面就會(huì)受到?jīng)_刷。隨著植被的減少,河流的沖刷會(huì)加劇。隨著洪水的沖刷,潰口的深度會(huì)越來(lái)越大,當(dāng)達(dá)到一定的穩(wěn)定程度后,潰口就會(huì)失去穩(wěn)定,然后慢慢向兩邊蔓延,最后形成一個(gè)梯形。

對(duì)于管涌潰口來(lái)說(shuō),在滲流過(guò)程中,水流作用會(huì)使土壤中的細(xì)顆粒被沖刷出粗大孔洞,在土壤中形成一個(gè)貫穿的滲流通道。隨著時(shí)間的推移,逐漸發(fā)展為暗渠式的潰口,當(dāng)侵蝕加劇時(shí),就會(huì)發(fā)生潰壩。此時(shí),潰口流量的計(jì)算公式如下:

Wa=M[2λ+(I-Iρ)/1+fL/E]0.5

(2)

式中:Wa為滲透通道的實(shí)際流量;M為潰口的橫斷面積;λ為常數(shù);I-Iρ為潰口的水頭;f為摩擦因子;L為滲透管道的實(shí)際長(zhǎng)度;E為滲透管道的實(shí)際直徑。

確定潰決的相關(guān)模式后,可以選取BREACH模型運(yùn)算所需參數(shù)。研究選取的模型運(yùn)算參數(shù)見(jiàn)圖2。

圖2 BREACH模型運(yùn)算所需參數(shù)示意圖

由圖2可知,BREACH模型運(yùn)算時(shí),所需參數(shù)包含壩頂與壩底的高度、壩高的最大值、溢洪道的頂部高度、壩頂?shù)膶?shí)際長(zhǎng)度與寬度、上游與下游壩坡比、壩體材料的粒徑和黏聚力、校核洪峰流量和干流的平均坡降、壩體材料的濕密度和摩擦角。

1.2 MIKE21模型原理及構(gòu)建分析

在河流流量計(jì)算的基礎(chǔ)上,研究選擇MIKE21軟件作為后續(xù)洪水演進(jìn)的數(shù)值模擬軟件。MIKE21是一種分析河流、河口等的二維化仿真模擬工具[13-14],其特點(diǎn)見(jiàn)圖3。

圖3 MIKE21軟件技術(shù)的特點(diǎn)示意圖

由圖3可知,MIKE21軟件具備無(wú)障礙性,即具備完整的人機(jī)交互功能,容易使用;高效性,即擁有快速可靠的模擬引擎;便捷性,即具備相應(yīng)的支持軟件用于數(shù)據(jù)處理和分析;廣泛性,即具備滿(mǎn)足各種河口模擬需求的相應(yīng)模塊;兼容性,即用于在第三方程序中處理模擬數(shù)據(jù)的地理信息系統(tǒng)集成和工具;靈活性,即多種計(jì)算網(wǎng)格,模塊和許可選項(xiàng)保證用戶(hù)可以根據(jù)實(shí)際需求選擇模型;公認(rèn)性,即具備長(zhǎng)達(dá)25年的相關(guān)記錄以及在全球廣泛應(yīng)用所證明的技術(shù)。

MIKE21的水動(dòng)力數(shù)學(xué)模型以二維非穩(wěn)定流動(dòng)為基本控制方程,其基本假定是以水為不可壓縮水體,并且壓力在水深上是按靜壓力分布的[15]。其中,二維非穩(wěn)定流動(dòng)基礎(chǔ)控制中的空間離散常用的是有限容積方法,即控制容積方法。在實(shí)際的二維模型中,由于模擬區(qū)的邊界往往是不規(guī)則的,因此其把二維計(jì)算區(qū)分成一些無(wú)組織的或不規(guī)則的小單元,再對(duì)各單元進(jìn)行水量和動(dòng)量的均衡,從而得出各計(jì)算單元的水位和流量。

二維模型一般采用結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。其中,非結(jié)構(gòu)化網(wǎng)格可以使網(wǎng)格合理分布,適合于復(fù)雜地形的網(wǎng)格劃分。研究針對(duì)河北省某縣的邊界條件,將其分為不規(guī)則的非結(jié)構(gòu)化網(wǎng)格。網(wǎng)格的大小隨著地形的變化而變化,對(duì)高程變化小、邊界規(guī)則的地區(qū)劃分范圍更廣;對(duì)高程變化大、邊界復(fù)雜的地區(qū)劃分的面積會(huì)相應(yīng)減少,也會(huì)進(jìn)行局部加密,以增加計(jì)算效率和準(zhǔn)確度。

因此,研究利用MIKE21軟件對(duì)河北省邢臺(tái)水庫(kù)的下游城區(qū)進(jìn)行二維模型計(jì)算。計(jì)算面積21.88km2,而網(wǎng)格剖分的相關(guān)部分選擇不規(guī)則的三角形網(wǎng)絡(luò),通過(guò)剖分得到25 000多個(gè)網(wǎng)孔,平均網(wǎng)孔面積870×10-6km2,最大網(wǎng)孔為5 260×10-6km2,最小為12×10-5km2。MIKE21模型一般將網(wǎng)孔尺寸控制在0.05km2以?xún)?nèi),因此研究劃分的網(wǎng)孔尺寸滿(mǎn)足邢臺(tái)水庫(kù)潰壩洪水模擬的精度要求。

在潰壩洪水模擬中,邊界條件對(duì)計(jì)算結(jié)果的精確度有較大影響。在二維模型計(jì)算中,存在著開(kāi)放邊界和封閉邊界兩種情況,前者通常為陸邊界,后者通常為水邊界。在實(shí)際的潰壩洪水?dāng)?shù)值模擬中,依據(jù)網(wǎng)格劃分水庫(kù)潰壩洪水演進(jìn)的數(shù)值模型流程見(jiàn)圖4。

圖4 網(wǎng)格劃分方法下潰壩洪水演進(jìn)模型流程示意圖

由圖4可知,構(gòu)建潰壩洪水演進(jìn)模型流程首先是對(duì)水庫(kù)下游相關(guān)區(qū)域地形進(jìn)行必要的處理,研究選擇的處理軟件為地理信息系統(tǒng)(Geographic Information System,GIS)。其次是利用MIKE21軟件內(nèi)部的相關(guān)建模工具,對(duì)網(wǎng)格進(jìn)行相應(yīng)的劃分;之后將GIS提出的相關(guān)高程點(diǎn)導(dǎo)入網(wǎng)格文件,進(jìn)行數(shù)字高程插值,并利用MIKE21生成模擬區(qū)域的實(shí)際地形圖;然后對(duì)BREACH和MIKE21模型的相關(guān)參數(shù)進(jìn)行設(shè)置;最后將BREACH的數(shù)值計(jì)算結(jié)果視為邊界條件,導(dǎo)入MIKE21模型中進(jìn)行耦合,以此完成整個(gè)模型的實(shí)際構(gòu)建。

2 耦合模型在潰壩洪水?dāng)?shù)值模擬中的應(yīng)用

為了驗(yàn)證BREACH與MIKE21的耦合模型在潰壩洪水?dāng)?shù)值模擬中的應(yīng)用效果,研究選取河北省邢臺(tái)水庫(kù)作為研究對(duì)象,將模型應(yīng)用在該水庫(kù)某年洪水事故下的潰壩洪水?dāng)?shù)值模擬中,并將該水庫(kù)下游水庫(kù)設(shè)定為Q。試驗(yàn)前,依據(jù)潰決原因,將潰決模式分為因洪水漫頂潰決的邢臺(tái)水庫(kù)和Q水庫(kù)連潰和Q水庫(kù)單潰(方案1和方案7)、因地震潰決(方案2和方案3)和管涌潰決(方案4至方案6)的邢臺(tái)水庫(kù)單潰。其中,因?yàn)闈Q水位不同又可以將邢臺(tái)水庫(kù)單潰分為5種方案,總計(jì)7種方案。因此,利用耦合模型模擬的不同方案下邢臺(tái)Q上下水庫(kù)潰口處流量過(guò)程曲線(xiàn)見(jiàn)圖5。

由圖5可知,不同方案下的潰口洪峰流量均大于10 000年一遇的入庫(kù)洪峰高達(dá)2.31×104m3/s的流量。其中,方案1最大的實(shí)際下泄流量為18.39×104m3/s;方案2為15.31×104m3/s;方案3為13.30×104m3/s;方案4為8.41×104m3/s;方案5為2.37×104m3/s;方案6為2.99×104m3/s;方案7為2.71×104m3/s。

圖5 不同方案下的邢臺(tái)和Q水庫(kù)潰口處流量過(guò)程曲線(xiàn)示意圖

綜合來(lái)看,不同方案下的潰壩洪水是以單一的波浪形式向下游擴(kuò)散的,在潰壩后,潰壩的下泄流量很快達(dá)到一個(gè)高峰,隨后逐漸下降。根據(jù)各種方案的對(duì)比發(fā)現(xiàn),潰壩潰決的水位越高,潰壩的洪峰流量越大;洪水決堤時(shí)間越短,洪水的沖刷流量越大。隨著決口水位的升高,泄洪時(shí)間也在逐漸延長(zhǎng)。在邢臺(tái)水庫(kù)和Q水庫(kù)連崩時(shí),Q水庫(kù)的洪峰流量依賴(lài)于邢臺(tái)水庫(kù)的潰決流量。邢臺(tái)水庫(kù)一旦出現(xiàn)單潰,豐寧上下水庫(kù)就會(huì)同時(shí)潰堤。

在此基礎(chǔ)上,針對(duì)各方案下邢臺(tái)水庫(kù)位置至Q水庫(kù)以及Q水庫(kù)下游沿線(xiàn)潰壩洪水的實(shí)際代表斷面進(jìn)行分析。由于篇幅所限,研究選取沿邢臺(tái)七里河的5個(gè)重要斷面地點(diǎn),分別為兩水庫(kù)之間的兩個(gè)地點(diǎn)以及Q水庫(kù)下游沿河3個(gè)地點(diǎn),分別用A、B、C、D、E表示。實(shí)際的洪峰流量變化趨勢(shì)見(jiàn)圖6。

由圖6可知,方案1決口洪峰向Q水庫(kù)的洪峰流量為1.46×104m3/s,比邢臺(tái)水庫(kù)的泄洪速度下降約20%;洪峰向地點(diǎn)C匯入點(diǎn)的洪水最大流量為1.01×104m3/s,比邢臺(tái)水庫(kù)壩址下降約44%。地點(diǎn)E的洪峰流量達(dá)到7×104m3/s,比邢臺(tái)水庫(kù)的洪水下降約57%。而地震潰決模式的方案2和方案3從邢臺(tái)水庫(kù)向Q水庫(kù)傳播時(shí),流量分別降低約18%和17.9%;流至地點(diǎn)E時(shí)分別降低約60%和59%。管涌潰決模式的方案4和方案5從邢臺(tái)水庫(kù)向Q水庫(kù)傳播時(shí),流量分別為7×104和2×104m3/s,分別降低約16%和4%;最后的方案6和方案7表現(xiàn)出同樣的效果。

圖6 各方案下潰壩洪水沿線(xiàn)洪峰流量的變化趨勢(shì)

綜合來(lái)看,在同一潰決水位條件下,潰壩洪峰流量隨長(zhǎng)度的增大而逐漸減小。而同一潰決時(shí)間段內(nèi),潰決時(shí)間隨潰決水位下降而減小;在同一潰決水位上,洪水消退時(shí)間愈長(zhǎng),洪水流向下游的速度愈慢。另外,潰壩洪水向下游演進(jìn)的實(shí)際過(guò)程中,不同方案下的沿線(xiàn)主要斷面最高水位都有所不同。研究利用耦合模型對(duì)其進(jìn)行數(shù)值模擬,見(jiàn)圖7。

圖7 不同方案下河流沿線(xiàn)最大洪水水位的變化趨勢(shì)示意圖

由圖7可知,方案1至方案7在Q水庫(kù)時(shí)的最高水位分別為306、300、297、286、272、270以及270m。到地點(diǎn)C時(shí)的最高水位分別為279、274、271、261、251、249以及247m。綜合來(lái)看,決口水位和決口時(shí)間是決定河道最高水位的重要因素,隨著決口水位的升高,河道的最高水位也隨之升高。在同一決口水位下,決口時(shí)間愈短,最高水位愈高。

研究利用耦合模型,在5個(gè)重要斷面中選取其中兩個(gè)斷面,并對(duì)其水位與時(shí)間的過(guò)程進(jìn)行分析,結(jié)果見(jiàn)圖8。

圖8 不同方案下兩個(gè)重要斷面的水位特征

由圖8可知,方案1的潰壩洪水水位始終最高,其中A地最高為306m,B地最高為290m。而方案7的潰壩洪水水位始終最低,其中A地最高為265m,B地最高為263m。綜合來(lái)看,各斷面洪水的水位隨著時(shí)間不斷向前推移而瞬間上漲,達(dá)到最高水位之后,又緩慢降低,直到穩(wěn)定水位。

3 結(jié) 論

為了在水庫(kù)遭遇潰壩之前制定合理的防洪避險(xiǎn)預(yù)案,本文將BREACH模型與MIKE21融合成耦合模型,并將其運(yùn)用在河北邢臺(tái)水庫(kù)潰壩的洪水?dāng)?shù)值模擬試驗(yàn)中。結(jié)果表明,不同方案下的潰口洪峰流量均大于高達(dá)2.31×104m3/s的流量;方案1決口洪峰向Q水庫(kù)的洪峰流量為1.46×104m3/s,比邢臺(tái)水庫(kù)的泄洪速度下降約20%;7個(gè)方案在Q水庫(kù)時(shí)的最高水位為方案1的306m;方案1的潰壩洪水水位始終最高。綜合來(lái)看,本文提出的耦合模型,在模擬河北省邢臺(tái)潰壩洪水的數(shù)值演變上具備較高的有效性。