MIKE軟件技術在水庫潰壩洪水數值模擬中的應用研究

管 靜

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

0 引 言

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

在此背景下,本文引入土壩失事洪水演算模型(BREACH)與二維數值工程軟件(MIKE21),并將二者組合構建成耦合模型,同時將其運用在河北邢臺水庫潰壩的洪水數值模擬應用中。其目的是為邢臺水庫在形成潰壩之前制定合理的防洪避險方案,并通過預測潰壩對下游的影響,為下游制定進一步的防洪計劃和轉移方案提供理論依據。此外,通過將耦合模型應用在潰壩洪水數值模擬過程中,制定7種不同的方案,既能全方位分析不同情況的影響,又具備一定的創新性。

1 基于MIKE技術的水庫潰壩洪水數值模擬分析

1.1 水庫潰口流量計算模型研究

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

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

BREACH模型可以細分為漫頂潰決和管涌潰決。對于漫頂潰決來講,當坡面沒有植被覆蓋時,由于水流的沖刷,在下游坡面上逐漸形成一個矩形的溝道。此時,河流水流流量計算公式如下:

Pa=3A0(h-hb)

(1)

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

因此,水壩潰口的實際發展過程見圖1。

圖1 水壩潰口的發展示意圖

由圖1可知,在植被覆蓋的情況下,如果坡面上的水流速度超過允許速度,下游的坡面就會受到沖刷。隨著植被的減少,河流的沖刷會加劇。隨著洪水的沖刷,潰口的深度會越來越大,當達到一定的穩定程度后,潰口就會失去穩定,然后慢慢向兩邊蔓延,最后形成一個梯形。

對于管涌潰口來說,在滲流過程中,水流作用會使土壤中的細顆粒被沖刷出粗大孔洞,在土壤中形成一個貫穿的滲流通道。隨著時間的推移,逐漸發展為暗渠式的潰口,當侵蝕加劇時,就會發生潰壩。此時,潰口流量的計算公式如下:

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

(2)

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

確定潰決的相關模式后,可以選取BREACH模型運算所需參數。研究選取的模型運算參數見圖2。

圖2 BREACH模型運算所需參數示意圖

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

1.2 MIKE21模型原理及構建分析

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

圖3 MIKE21軟件技術的特點示意圖

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

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

二維模型一般采用結構化網格和非結構化網格。其中,非結構化網格可以使網格合理分布,適合于復雜地形的網格劃分。研究針對河北省某縣的邊界條件,將其分為不規則的非結構化網格。網格的大小隨著地形的變化而變化,對高程變化小、邊界規則的地區劃分范圍更廣;對高程變化大、邊界復雜的地區劃分的面積會相應減少,也會進行局部加密,以增加計算效率和準確度。

因此,研究利用MIKE21軟件對河北省邢臺水庫的下游城區進行二維模型計算。計算面積21.88km2,而網格剖分的相關部分選擇不規則的三角形網絡,通過剖分得到25 000多個網孔,平均網孔面積870×10-6km2,最大網孔為5 260×10-6km2,最小為12×10-5km2。MIKE21模型一般將網孔尺寸控制在0.05km2以內,因此研究劃分的網孔尺寸滿足邢臺水庫潰壩洪水模擬的精度要求。

在潰壩洪水模擬中,邊界條件對計算結果的精確度有較大影響。在二維模型計算中,存在著開放邊界和封閉邊界兩種情況,前者通常為陸邊界,后者通常為水邊界。在實際的潰壩洪水數值模擬中,依據網格劃分水庫潰壩洪水演進的數值模型流程見圖4。

圖4 網格劃分方法下潰壩洪水演進模型流程示意圖

由圖4可知,構建潰壩洪水演進模型流程首先是對水庫下游相關區域地形進行必要的處理,研究選擇的處理軟件為地理信息系統(Geographic Information System,GIS)。其次是利用MIKE21軟件內部的相關建模工具,對網格進行相應的劃分;之后將GIS提出的相關高程點導入網格文件,進行數字高程插值,并利用MIKE21生成模擬區域的實際地形圖;然后對BREACH和MIKE21模型的相關參數進行設置;最后將BREACH的數值計算結果視為邊界條件,導入MIKE21模型中進行耦合,以此完成整個模型的實際構建。

2 耦合模型在潰壩洪水數值模擬中的應用

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

由圖5可知,不同方案下的潰口洪峰流量均大于10 000年一遇的入庫洪峰高達2.31×104m3/s的流量。其中,方案1最大的實際下泄流量為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 不同方案下的邢臺和Q水庫潰口處流量過程曲線示意圖

綜合來看,不同方案下的潰壩洪水是以單一的波浪形式向下游擴散的,在潰壩后,潰壩的下泄流量很快達到一個高峰,隨后逐漸下降。根據各種方案的對比發現,潰壩潰決的水位越高,潰壩的洪峰流量越大;洪水決堤時間越短,洪水的沖刷流量越大。隨著決口水位的升高,泄洪時間也在逐漸延長。在邢臺水庫和Q水庫連崩時,Q水庫的洪峰流量依賴于邢臺水庫的潰決流量。邢臺水庫一旦出現單潰,豐寧上下水庫就會同時潰堤。

在此基礎上,針對各方案下邢臺水庫位置至Q水庫以及Q水庫下游沿線潰壩洪水的實際代表斷面進行分析。由于篇幅所限,研究選取沿邢臺七里河的5個重要斷面地點,分別為兩水庫之間的兩個地點以及Q水庫下游沿河3個地點,分別用A、B、C、D、E表示。實際的洪峰流量變化趨勢見圖6。

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

圖6 各方案下潰壩洪水沿線洪峰流量的變化趨勢

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

圖7 不同方案下河流沿線最大洪水水位的變化趨勢示意圖

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

研究利用耦合模型,在5個重要斷面中選取其中兩個斷面,并對其水位與時間的過程進行分析,結果見圖8。

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

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

3 結 論

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