吉林省飲馬河洪水演進數(shù)值模擬研究

李 航，丁 曼，徐貴新

(吉林省水利水電勘測設計研究院，長春 130021)

飲馬河是吉林省第二松花江左岸最大支流，發(fā)源于四平市伊通縣地局子鄉(xiāng)老爺嶺東南，流經(jīng)伊通縣、磐石市、永吉市、雙陽區(qū)、九臺區(qū)、德惠市、農(nóng)安縣和長春市二道區(qū)等8個市縣(區(qū))，于農(nóng)安縣靠山鎮(zhèn)紅石壘屯匯入第二松花江，集水面積18 247 km2，其中石頭口門水庫以上集水面積4 944 km2。飲馬河河道長度387.5 km，流域平均高程300 m左右。伊通河流域其上游屬丘陵崗地，下游地形波狀起伏，兩岸河谷較開闊，地勢較高部分大都開墾為旱田，河谷平地大部分開墾為水田。流域水系發(fā)育、支流較多，較大支流依次有伊通河、雙陽河、霧開河、岔路河、三道溝河和小南河等[1]。

分析范圍為飲馬河石頭口門水庫至飲馬河河口，包括沿岸的防洪保護區(qū)。其中，一維河道模型的構建范圍為石頭口門水庫到飲馬河河口，范圍內(nèi)的主干河道長度為110 km；二維水流模型的構建范圍為飲馬河沿岸的防洪保護區(qū)。

圖1 飲馬河(石頭口門水庫以下)水系示意圖

1 洪水分析方法及模型構建

1.1 洪水分析方法

飲馬河石頭口門水庫放流與區(qū)間小南河、霧開河與伊通河等支流組成編制區(qū)域的洪水來源，干流來水控制站為石頭口門水庫站和德惠站，支流來水控制站為九臺站、農(nóng)安站等。

采用MIKE ZERO系列洪水模擬水動力軟件進行洪水分析，包含河道洪水演進、潰堤洪水以及淹沒區(qū)洪水演進，分別建立MIKE11一維非恒定流模型模擬飲馬河河道洪水演進過程，DAMBRK方法計算堤防潰堤洪水，MIKE21二維非恒定流模型模擬飲馬河防洪保護區(qū)的洪水演進過程[2]。由于二維水流模擬的上邊界為堤防潰口處的流量過程，而潰口流量又取決于潰口上下游(河道與淹沒區(qū))的水位差，因此將上述3種模型相互耦合、聯(lián)立求解，系統(tǒng)模擬潰堤后對防洪保護區(qū)的影響。

1.2 模型原理

1.2.1 一維水動力學模型原理

圣維南(Saint-Venant)方程組為MIKE11一維非恒定流模型中的河道洪水運動的基本方程：

動量方程：

式中：Q為斷面流量，m3/s；A為過水斷面面積，m2；q為源匯的單寬流量，m2/s；x為距離坐標，m；t為時間坐標，s；h為水位，m；C為謝才系數(shù)；R為水力半徑，m；g為重力加速度，m/s2；α為動量校正系數(shù)。

1.2.2MIKE21二維水動力模型

二維非恒定流方程如下，為平面二維水流模型，由連續(xù)方程及動量方程組成：

式中：h為水深，m；Z為水位，m；u為X方向沿垂線平均的水平流速分量，m/s；v為Y方向沿垂線平均的水平流速分量，m/s；g為重力加速度，m/s2；n為糙率；q為源匯項，m/s。

1.3 模型構建

1.3.1 網(wǎng)格剖分

收集編制范圍1∶10 000矢量地形圖及數(shù)字高程，確定模型中防洪保護區(qū)分析范圍，并設置開邊界和閉邊界，開邊界為水流進出邊界，閉邊界為阻水邊界?？紤]無結構網(wǎng)格適應能力強并且靈活的特點，為提高邊界處模擬精度，采用無結構三角形網(wǎng)格剖分[3]，控制三角形網(wǎng)格最大面積不超過0.01 km2，且網(wǎng)格大小均勻，盡量剖分成等邊三角形。見圖2、圖3。

圖2 網(wǎng)格劃分示意圖

圖3 局部網(wǎng)格劃分示意圖

1.3.2 計算時間和步長

河道斷面間距、連接位置以及參數(shù)驗證的水文站點位置等因素與一維河道模型空間步長設定息息相關，飲馬河斷面間距最小為50 m，作為MIKE模型中水位計算點間距，為了保證模型穩(wěn)定準確運行，此次時間步長設定為30 s。

1.3.3 邊界條件設置

根據(jù)編制范圍河網(wǎng)分布，上邊界位于石頭口門壩址，根據(jù)石頭口門水庫調(diào)度規(guī)則，調(diào)節(jié)天然入庫洪水過程線得到水庫泄流過程，下邊界位于飲馬河河口，根據(jù)河口處實測斷面分析計算得到水位-流量關系。區(qū)間小南河、霧開河與伊通河等支流作為點源以洪水過程線形式加載。

1.3.4 模型參數(shù)選取與驗證

根據(jù)飲馬河各斷面河道河床的實際組成情況，并通過石頭口門站、德惠站洪痕調(diào)查以及H-Q關系曲線反推糙率驗證而得到的糙率成果，主槽值選用0.03，邊灘選用0.05。

將防洪保護區(qū)按照土地利用圖分為6種地表類型，各類下墊面糙率參考《洪水風險圖編制導則》有關數(shù)值，并根據(jù)典型區(qū)域?qū)崪y大洪水淹沒資料進行分析率定。見表1。

表1 飲馬河防洪保護區(qū)糙率表

1.3.5 模型驗證

以飲馬河石頭口門水庫2010年7月21日實測泄流為上邊界，結合水文年鑒中洪水要素摘要實測流量，模擬大汛期洪水，進行模型參數(shù)驗證。通過調(diào)整飲馬河各斷面的河床主槽及邊灘的糙率參數(shù)，使干流區(qū)間內(nèi)的德惠水文站處經(jīng)模型模擬得到的流量和水位與實測值接近吻合。當洪峰流量Qm和水位H的誤差分別控制在10%和20 cm以內(nèi)[4]，說明河道參數(shù)取值較為合理。見圖4。

利用已經(jīng)參數(shù)驗證后的一維河道模型，采用飲馬河干流及各支流50年一遇設計洪水(以1953年作為典型年，干流洪峰流量Qm為700 m3/s)作為邊界條件，分析計算飲馬河石頭口門水庫壩下至河口段河道水面線，并與已通過上級水利主管部門審查的50年一遇設計水面線成果進行對比(圖5)，模擬水面線成果與堤防設計成果平均誤差為10 cm，最大誤差為21.4 cm，驗證結果符合要求。

圖4 德惠站模擬值與實測值對比圖(左圖為流量，右圖為水位)

圖5 模擬與設計水面線對比圖(P=2%)

2 洪水分析成果

2.1洪水分析計算

以陳家崴子潰口20年一遇洪水方案為例，飲馬河干流按20年一遇洪水控制，各支流相應洪水組合，洪水過程2015年8月2日0時至9月20日6時。8月28日8時，德惠市農(nóng)村段陳家崴子險工處河道達到特征水位，堤防開始發(fā)生潰決，最大潰口寬度50 m，潰口最大流量Qm=297 m3/s。8月24日12時，飲馬河下游里程75 606 m處有小部分水漫出，其余堤防均未發(fā)生漫堤現(xiàn)象。至洪水過程結束，共進入飲馬河防洪保護區(qū)總水量4 451×104m3，最大淹沒范圍61.22 km2，最大淹沒水深達4.42 m，最大洪水流速達1.29 m/s。見圖6、圖7。

2.2 水量平衡分析

2.2.1 淹沒區(qū)水量平衡分析

計算模型中的進水量、出流量以及防洪保護區(qū)內(nèi)淹沒的水量，并將入流量減去出流量后得到的差值，與防洪保護區(qū)內(nèi)淹沒的水量對比，若誤差在1%以內(nèi)，說明洪水模擬結果合理。經(jīng)計算，陳家崴子潰口方案進水量11 209.78×104m3，出水量6 759.26×104m3，最終水量4 451.41×104m3，相對誤差0.02%。

圖6 飲馬河陳家崴子潰口流量過程圖

圖7 飲馬河陳家崴子潰口最大淹沒水深示意圖

2.2.2 耦合模型水量平衡分析

在一、二維耦合模型計算中，潰堤模型作為橋梁，將一維水動力模型水量傳遞給二維水動力模型。由于是三者耦合計算，如模型設置不合理，會造成水量缺失現(xiàn)象，影響結果的可信度。通過對比一維模型輸入二維模型水量與二維模型接收一維模型水量情況檢驗計算合理性，陳家崴子潰口發(fā)生20年一遇洪水時，一維水動力模型河道出流水量11 209.78×104m3，二維水動力模型邊界入流水量11 209.78×104m3，水量相差為0，說明模型計算合理[5]。

3 結 語

通過MIKE系列軟件構建一、二維水動力模型，可以直觀展示洪水在發(fā)生潰堤后產(chǎn)生的最大淹沒范圍、最大淹沒水深及最大洪水流速等信息，同時也可展示洪水演進動態(tài)過程，能夠提取洪水到達時間、洪水淹沒歷時等洪水風險要素。本文以飲馬河陳家崴子潰口20年一遇洪水方案為例，介紹了MIKE11和MIKE21模型原理，模型構建中的網(wǎng)格剖分、設置邊界條件、模型參數(shù)驗證等，建立具有較高模擬精度的水動力耦合模型，模擬成果經(jīng)水量平衡分析后合理可靠，成果可為后續(xù)的洪災損失評估和防洪減災提供技術支撐。