基于MIKE21 模型的黃河四排口河段洪水演進數值模擬

鎖曉南 李春光＊，2 尚彥祥 王 悅

（1.寧夏大學土木與水利工程學院，寧夏 銀川 750021；2.北方民族大學土木工程學院，寧夏 銀川 750021）

0 引言

目前，對于河道洪水演進的數學模型主要是基于一維、二維水動力數學模型。 王揚等[1]基于一維、二維水動力耦合模型建立了韓江南北堤防洪保護區潰壩洪水演進數學模型。 周躍華等[2]基于GIS 與一維水動力模型對陶樂防洪保護區漫溢洪水進行了風險分析。孫玲玲[3]通過建立二維非恒定流數學模型模擬了黃壁莊水庫洪水演進過程。 戚藍等[4]基于高精度DEM 數據建立三亞市主城區暴雨、洪水與風暴潮多元耦合的精細化洪澇分析模型。 王俊琿等[5]通過建立高分辨率數值模型與三維可視化耦合技術對陜西灃西新城進行了洪澇過程分析。目前來說一維水動力數學模型計算簡單方便，模擬精度較高，在水庫洪水演進模擬中應用較多，但對于復雜的游蕩性河段來說其計算可靠性較差；基于二維水動力數學模型目前不僅應用于洪水演進中， 在河床變形與泥沙運移中也應用較廣，對復雜河道的計算精度高、可靠性較準確，能夠有效的模擬復雜動邊界條件下的洪水演進過程。

本文主要基于二維水動力數學模型， 采用有限體積法對黃河四排口河段進行洪水演進的模擬計算，分析大流量下河道的水動力場變化狀況，以及河道的沖刷和河床變形狀況， 對今后該河道的治理提供理論依據。

1 研究區域概況

黃河四排口河段位于寧夏石嘴山市平羅縣境內，該段屬于平原游蕩非穩定分叉型河道，河寬1 800～6 000 m，平均寬為3 300 m，主槽寬500～-1 000 m，平均寬為650 m,河道比降為0.18‰,河床糙率為0.016～0.027，彎曲率為1.23。 平面上呈寬窄相間的藕節狀分布，河道寬淺，水流散亂，沙洲密布，河床組成為細沙、粉細砂，河床抗沖性差，沖淤變化較大，主流游蕩擺動較為劇烈，兩岸主流頂沖點不定，經常出現險情。其左岸為銀川平原的組成部分，土地肥沃，是引黃老灌區，耕地面積30 萬余畝，右岸屬鄂爾多斯平原，沙質土壤，耕地面積較少。 由于該段是多年的老險工河段，歷史上每次遇到大洪水均發生頂沖榻岸、險情不斷，搶險投入較大，嚴重威脅平羅灌區。 為了從根本上消除洪水威脅，歸順黃河流路，寧夏自治區政府于2017—2018年對黃河四排口河段進行了‘裁彎曲直’整治工程，此次整治工程主要是改變四排口河段的水流方向，為此修建了26 座丁壩， 并構筑了聯壩將所有丁壩連接起來，且開挖了一條長約1.7 km，直線段上開口82 m，下開口40 m，深約5 m 的引河，以改變水流方向。

本文主要研究‘裁彎曲直’工程后大流量洪水對該段河道的沖淤狀況，針對此河段本文主要采用二維水動力MIKE 21 模型， 基于2018 年測量的地理數據構建黃河四排口河段水動力模型，選取2018 年10 月21 日與11 月25 日所測量的數據作為數據來源，對所建立的水動力數學模型進行驗證與率定，設置兩種工況研究該河道洪水期水流對河床的沖刷狀況，為今后該河段的治理提供建設性的意見。

2 四排口二維水動力數學模型的構建

2.1 數學模型的建立

本文所研究的黃河四排口河段的水平尺度遠遠大于垂向尺度[6]，重力是水流運動的主要驅動力，流速等參數沿水平方向上的變化遠遠大于沿垂直方向上的變化，水壓力接近靜壓分布等。 所以基于以上條件下建立該河段的沿水深方向的平面二維水動力數學模型。

連續方程：

其中：z 為水位；h 為水深；u 和v 分別為x 和y 方向上的流速；Vt為水流紊動粘性系數；g 為重力加速度；n 為曼寧系數；ρ 為水流的密度；τsx和τsy分別為x和y 方向的水面風應力；f0為柯氏力系數,柯氏力項一般可以忽略；v 為水流運動粘性系數；k 為紊動動能；ε為紊動動能耗散率；SK和Sε分別為k 方程和ε 方程中的源項。

2.2 地形建立及網格劃分

本文采用2018 年10 月21 日實測地形資料進行地形插值，其地形圖如下圖1 所示，采用非結構三角形網格對地形圖進行剖分，所允許剖分的非結構三角形網格最大面積不超過300 m2, 共計得到7 401 個網格節點，生成14 042 個平面網格單元，生成的網格模型圖見下圖2、圖3 所示。

圖1 四排口地形圖

2.3 模型初始條件和邊界條件的設置

2.3.1 時間步長

該模型采用的時間步長為180 s， 時間步數為16 800 步。

2.3.2 CFL 數

為了協調不同大小的網格，使其CFL 數滿足該模型的要求（CFL<1），對該模型使用變時間步長0.01～120 s。

圖2 四排口河段計算網格圖

圖3 丁壩段網格剖分圖

2.3.3 床面阻力系數

本文主要采用常見的對模型進行率定的方法，得到該研究區域的曼寧系數，因此本文的模型取值為32。

2.3.4 初始條件

初始水位為零，流速也為零；

2.3.5 邊界條件

進口邊界條件，黃河四排口河段上游主要設置為流量邊界；出口邊界條件，河段下游主要設置為水位邊界條件；河床邊界主要采用無滑動的陸地邊界。

2.4 水動力模型的驗證

本文采用2018 年10 月21 日現場實測的四排口河段水文數據（進口流量：948.62 m3/s、出口水位：1 097.35 m、沿程水位為h）對該河道所建立的模型進行率定與驗證，為了與計算結果方便對比，在該研究區域內選擇了13 個典型斷面作為監測點， 研究斷面的實測水位、模擬水位、以及水位差，對模擬結果進行驗證。 通過比對模擬水位與實測水位得出水位差，該模型下模擬出的各段面模擬水位與實測水位大致接近，計算得出其水位差的范圍為-0.016～0.041 m，相對誤差均控制在2%以內。該誤差在計算可接受范圍內，因此模型可用于黃河四排口河段洪水演進模擬計算。

3 洪水演進模擬結果分析

在水動力模型驗證計算的基礎上， 本文主要開展了兩組工況（平水期與豐水期）的計算，兩種工況的來流過程是一致的。 根據水文資料統計，近五十年來， 寧夏出現代表性的洪水年份主要為2012 年和1981 年。 2012 年洪水，青銅峽水文站測得最大洪峰流量為3 070 m3/s；1981 年洪水， 青銅峽水文站測得最大洪峰流量為6 040 m3/s。 因此本文豐水期的流量采用1981 年測得的最大洪峰流量6 040 m3/s 為計算流量，平水期的流量采用2018 年7 月17 日實測流量1 447 m3/s 為計算流量，出口斷面水位以2018 年7 月17 日實測水位為基準，每級流量均歷時30 d。 該模擬結果主要關注了河段的水動力場，局部河段的變形狀況。

該河段由于‘裁彎曲直’工程后河寬的大小沿水流運動的變化較大， 原始河道河寬平均約為1.43 km，新開挖引河段河寬約為135 m，直線段上開口為82 m，下開口為40 m，河段在丁壩末端急劇收縮，使得過流斷面減少，流速增大，形成了一種“瓶頸狀”的特殊河段。計算結果表明平水期工況下， 水流的運動較為平緩，過度河段較能夠過流；豐水期（洪水）工況下，由于河段上游左岸有一系列丁壩構成的聯壩，造成該河段的局部阻力增大，使得泄流能力下降，引發了雍水現象，并且上游來水在丁壩的挑流作用下，將上游來水沿丁壩迎水面壩坡挑至河道中央，致使下泄主流方向明顯向右岸發生偏轉，并在單個丁壩的背水面出現了明顯的壩后回流現象，但由于丁壩背水面水位較低且此處流速較小因此對此處影響較小，見下圖5。

圖4 河段流場分布圖

在豐水期（洪水）工況下，圖5 中可以明顯看出丁壩處的水流速度很大且伴隨著長時間的沖刷，會嚴重影響丁壩的安全性與穩定性；圖6 中由于上游一系列的丁壩使得下游左岸出現了一個長半徑約280 m，寬半徑約80 m 的逆時針運動相對規則的橢圓形漩渦，相比于平水期的小漩渦， 該工況下的漩渦明顯增大，其形成原因主要為流量增大、水位上升、岸邊嚴重被沖刷，并且由于該段右岸河道邊坡較為陡峭，左岸河道比較順直，使得主流匯入該段先偏向右岸后在‘淹沒丁壩’ 的作用下流向河床高程更低的左岸區域，在這樣的流速情況下， 該段左岸必然會出現河床沖刷、邊岸崩塌后退現象，而右岸流速相對較低，在局部地區可能會發生大量淤積現象。

圖5 丁壩段流場分布圖

圖6 丁壩末端計算流場分布

4 結束語

本文主要建立了二維水動力數學模型并利用2018 年實測數據進行了驗證， 利用所建立的二維水動力數學模型對黃河四排口河段進行了洪水演進和河床變形的數值模擬，結果表明。

（1）該河段沿程水位高程模擬值與實測值較為接近，其水位差范圍為-0.016～0.041 m，相對誤差均控制在2%以內， 因此所建立的二維水動力數學模型適用于該河段的洪水演進模擬計算。

（2）洪水期工況下，水流對丁壩的沖刷劇烈，嚴重影響了丁壩的穩定性與安全性；且在丁壩末端左岸形成了一個長半徑約274.76 m，寬半徑約78.38 m 的逆時針運動相對規則的橢圓形漩渦，加劇了對該段的沖刷，可能會使得引河段水流向左岸遷移，導致引河段左岸發生大幅度的崩塌后退狀況。