渝河(隆德段)洪水數(shù)值模擬研究

楊東旭, 李金燕

(1.寧夏回族自治區(qū)水旱災(zāi)害防御中心,寧夏 銀川 750001;2.寧夏大學 土木與水利工程學院,寧夏 銀川 750021)

洪水分析模擬是指通過洪水計算軟件計算特定頻率洪水的淹沒風險特征信息，確定特定洪水重現(xiàn)期或特定方案下的洪水位、可能淹沒范圍和淹沒水深等；洪災(zāi)損失評估是指通過災(zāi)情損失評估模型，結(jié)合區(qū)域洪水分析模擬結(jié)果，據(jù)此計算特定方案的洪災(zāi)損失.

一維(1D)水動力學模型主要是對水流運動過程進行模擬，它可以將復雜河網(wǎng)的水位變化及流量變化過程快速、準確地模擬出來，但模擬的流動路線難以避免“一維線性”特征的出現(xiàn).而對于水平漫流過程[1]的模擬，主要通過二維(2D)水動力學模型來實現(xiàn)，它是將時間進行有序的間隔設(shè)定，對淹沒深度進行模擬計算.對于二維水動力學模型來說，它具有計算時間長，對地形數(shù)據(jù)的分辨率要求較高[2]等缺點.由此可見，一維、二維水動力學模型描述水流問題的側(cè)重點是不同的.一般來說，一維水動力學模型是用來模擬河道洪水演進的流動過程，二維水動力學模型則是模擬解決防洪保護區(qū)的洪水流動問題.本文將兩種模型進行側(cè)向耦合求解，不僅可以將潰堤、漫堤洪水的淹沒特征通過二維平面凸顯出來，還可以提高水動力學模型的計算效率和模擬精度[3]，有效減少網(wǎng)格數(shù)量，增加模型可靠性[4].對于分洪對河道水位產(chǎn)生的影響，以及河道水位由于某些因素下降對蓄滯洪區(qū)或者城市分洪量造成的影響等問題，該耦合模型系統(tǒng)內(nèi)部均可進行自動計算分析.

本文應(yīng)用洪水仿真、GIS、數(shù)據(jù)庫等技術(shù)，獲得防洪保護區(qū)內(nèi)洪水演進的可能途徑和洪水所造成的淹沒情況，便于防汛部門制定防洪調(diào)度決策，及時部署指揮防洪搶險救災(zāi)任務(wù)，補充和完善洪澇災(zāi)情評估，為災(zāi)后重建和土地開發(fā)利用提供管理依據(jù).

1 數(shù)值模型

MIKE11模塊可以很好地模擬分析一維河道洪水的演進情況，該模型計算穩(wěn)定且計算精度高，具有較強的可靠性，適用于河道水文要素的計算，例如河道任意時刻、任意斷面的水位觀測和流量計算.該模型還可用來模擬水庫調(diào)度、蓄滯洪區(qū)運用，分析潰口等構(gòu)筑物對河道洪水演進的影響.MIKE21能夠較為精準地計算該保護區(qū)內(nèi)的二維洪水演進情況，可以具體模擬出水流水平方向和垂直方向的變化.以下為模型的計算原理.

1.1 一維水動力學模型

一維水動力學模型采用MIKE11模塊軟件模擬.模擬時需構(gòu)建一維非恒定流Saint-Venant方程組，即水量方程和動量守恒方程[5]，通過Abbott六點隱式格式離散方程組，按序交替計算每個網(wǎng)格節(jié)點處的水位和流量.這種離散格式使得大步長計算時能在非常大的常數(shù)下保持穩(wěn)定，同時節(jié)省計算時間.

(1)

(2)

式中：B為平均寬度；Z為水位；t為時間；Q為流量；s為距離；q為旁側(cè)入流單寬流量；A為過水斷面面積；g為重力加速度；C為謝才系數(shù)；R為水力半徑；i為坡度.

1.2 二維水動力模型

二維水動力學模型采用MIKE21模塊軟件模擬，該模型的控制方程采用改進形式的二維淺水方程[6]

(3)

動量方程為

(4)

(5)

式中：H為水深；Z為水位；q為旁側(cè)入流單寬流量；t為時間；M與N分別為x和y方向的垂向平均單寬流量；u和v分別為垂向平均流速在x與y方向的分量；n為曼寧糙度系數(shù).

1.3 一維二維耦合模型(MIKE FLOOD)

2 渝河(隆德段)洪水模擬

圖1 渝河計算區(qū)建模范圍示意圖

2.1 河道斷面設(shè)置與網(wǎng)格剖分

2.1.1河道斷面設(shè)置 河道斷面是一維水動力學模型的重要基礎(chǔ)數(shù)據(jù)[10]，將已收集的河道斷面數(shù)據(jù)按照MIKE11的斷面輸入要求整編處理，然后通過內(nèi)插加密，提高計算精確度，使計算結(jié)果最大限度貼近實際值.對于河道形態(tài)變化顯著的地方以及涉及工程(橋、閘、壩、堰等)的位置，斷面需進一步加密處理.經(jīng)一維水動力學模型計算，渝河河段長度為23.65 km，共設(shè)置了265個河道斷面，河道布置見圖2.

圖2 渝河河道斷面布置

2.1.2網(wǎng)格剖分 將二維非恒定流方程作為網(wǎng)格剖分時的基本控制方程，對模擬范圍進行無結(jié)構(gòu)不規(guī)則網(wǎng)格劃分，要求最大的不規(guī)則網(wǎng)格面積不能大于0.5 km2.網(wǎng)格剖分時需要根據(jù)地形、歷史洪水淹沒范圍、防洪工程分布情況等因素劃定風險區(qū)域，以便于后期分析計算，因此,地形變化較大的區(qū)域以及重要地區(qū)的網(wǎng)格要適當加密.本文利用網(wǎng)格剖分將渝河區(qū)域的二維平面劃分成17萬個小網(wǎng)格，其中面積最小為20 m2，最大為800 m2，總的計算面積達到43.21 km2，網(wǎng)格剖分示意圖見圖3.

圖3 二維區(qū)域網(wǎng)格剖分示意圖

2.2 參數(shù)確定

2.2.1河道一維水動力模型計算參數(shù) 本文將河道粗糙率作為影響渝河河道一維水動力學模型分析精度的計算參數(shù)[11]，參考相關(guān)設(shè)計報告[12—13]，設(shè)定渝河河道綜合糙率為0.03.將16座過水路面作為河道建筑物處理，概化為堰，鑒于路面具有阻水作用，根據(jù)設(shè)計參數(shù)設(shè)定過水尺寸時不可忽略路面阻水因素.綜合考慮模型運算效率及模型本身的穩(wěn)定性等多種因素[14]，將該模型計算迭代步長設(shè)定為1 s.根據(jù)渝河歷史洪災(zāi)數(shù)據(jù)記載，結(jié)合渝河的堤防現(xiàn)狀，在渝河(三里店水庫至小河口段)設(shè)置2處潰口，將潰口處河道水位與堤防設(shè)計水位持平時設(shè)為堤防潰決時刻臨界點，即此時堤防瞬間潰決，潰口寬度為60 m.和平村潰口底高程為1 901.7 m，神林村潰口底高程為1 818.3 m.

2.2.2計算區(qū)二維水動力模型計算參數(shù) 參考相關(guān)文獻，將防洪保護區(qū)綜合糙率值設(shè)為0.06[15—16]，綜合考慮模型運算效率及模型本身的穩(wěn)定性等多種因素[17]，將該模型計算迭代步長設(shè)定為10 s，最小迭代步長設(shè)為0.01 s.

2.3 邊界條件

根據(jù)支流匯入情況，將渝河分為3個河段分別進行洪水風險分析，本文選擇的洪水分析量級為20 a一遇、50 a一遇和100 a一遇，在渝河(三里店水庫至小河口段)設(shè)置2處潰口.根據(jù)實際調(diào)研結(jié)果，渝河洪水潰堤之前，保護區(qū)內(nèi)初始水深為0 m.因此，根據(jù)干水深和濕水深理論，可將計算區(qū)域的網(wǎng)格設(shè)定為干單元.

2.3.2渝河計算區(qū)二維水動力學模型邊界條件 渝河計算區(qū)二維水動力模型外邊界具有閉邊界和開邊界兩種.閉邊界即為與外界無水流交換的邊界，開邊界為一維模型和二維模型水流動態(tài)交互的邊界，河道潰口側(cè)向建筑物耦合的二維區(qū)域邊界可看作模型開邊界.以地面為原點，將高于地面0.5 m的阻水建筑物或過水建筑物作為水動力學模型的內(nèi)部邊界，例如堤防、線狀地物、涵洞、橋梁等.一般來說，洪水演進是指可以正常通過涵洞、橋梁等過水建筑物或者具有沿程缺口的線狀建筑物，而對于堤防等阻水建筑物來說，洪水的演進是以漫溢形式通過的.河道一維模型將為計算區(qū)平面二維模型實時提供漫溢或潰口分流流量，作為固定時間步長內(nèi)二維模型的入流條件.

3 模擬結(jié)果與合理性分析

通過對渝河(隆德段)10 a一遇、20 a一遇、50 a一遇以及100 a一遇的洪水進行模擬，最終建立渝河(隆德段)的洪水預案庫.本文將以100 a一遇洪水為代表進行模擬結(jié)果分析.

3.1 模擬結(jié)果分析

3.1.2潰口分流分析 和平村、神林村潰口(上游臨近斷面)河道流量變化過程見圖4，潰口處河道水位變化過程見圖5，潰口分流流量過程見圖6.

圖4 渝河100 a一遇洪水在和平村、神林村潰口(上游臨近斷面)河道流量變化過程

圖5 渝河100 a一遇洪水在和平村、神林村潰口處河道水位變化過程

本文將河道水位與設(shè)計水位相等時設(shè)定為河道潰堤臨界時刻[15].由圖4和圖5分析可知，和平村潰口處起潰時刻河道流量為120.3 m3/s，河道洪水水位為1 901.9 m，分流歷時為5.75 h；神林村潰口處起潰時刻河道流量為141.7 m3/s，此時河道洪水水位為1 818.7 m，分流歷時為5.5 h.由圖6分析可知，和平村潰口起潰時刻分流流量為16.7 m3/s，此時分流洪峰流量達到63.1 m3/s，而神林村潰口由于臨近上游斷面，起潰時刻分流流量是和平村的2.1倍，約為35 m3/s，同時其分流洪峰流量也達到了118.9 m3/s.

圖6 渝河100 a一遇洪水在和平村、神林村潰口分流流量過程

3.1.3漫溢分流分析 通過對渝河河道百年一遇洪水進行一維模擬發(fā)現(xiàn)，共有18處發(fā)生洪水漫溢現(xiàn)象，其中16處為過水路面.漫溢上段漫溢位置6處，下段漫溢位置12處.

3.1.4防洪保護區(qū)洪水淹沒分析 洪水進入二維平面計算區(qū)域的形式有潰流和漫溢兩種形式，選取1, 2, 4, 6.5 h共4個時段分別對兩個村潰口洪水淹沒水深的分布情況(圖7)進行比較分析.

由圖7可知，洪水演進1 h，和平村潰口已起潰，洪水由潰口向保護區(qū)演進.由于過水路面使水位壅高，渝河部分過水路面處出現(xiàn)洪水漫溢現(xiàn)象，計算區(qū)內(nèi)淹沒面積1.13 km2.洪水演進2 h，神林村潰口已起潰，洪水由潰口向保護區(qū)演進.渝河上段漫溢洪水水深大部分在0.5 m以下，下段河道洪水繼續(xù)向外漫溢，多處水深在0.5～1 m之間，由于潰堤洪水疊加漫溢洪水，神林村附近水深達到1 m以上，洪水向南北兩側(cè)演進過程中受到國道312的阻擋，路前最大積水深1.2 m，計算區(qū)淹沒面積4.9 km2.洪水演進4 h，漫溢洪水繼續(xù)向下游演進，國道312阻水效果明顯，淹沒面積6.16 km2.由于河道來洪減少，河道內(nèi)水位降低，二維保護區(qū)下游洪水通過堤防較低處流向河道.最終，淹沒面積為6.22 km2，積水量300.67萬m3.受大洪水沖刷影響，沿渝河河道橋梁、護坡護岸及堤防工程、過水路面、引退水工程等工程設(shè)施將遭受嚴重損害.

圖7 洪水演進淹沒水深分布

3.2 合理性分析

3.2.1水量平衡分析 本文的水量平衡關(guān)系可通過渝河防洪保護區(qū)進洪量和區(qū)內(nèi)蓄水量是否相等得以驗證.針對渝河防洪保護區(qū)設(shè)定的6個洪水方案分別進行進洪量和蓄水量的對比分析(見表1)，結(jié)果表明渝河不同洪水方案下進洪量等于蓄水量，誤差為0 m3，滿足水量平衡要求.

表1 渝河防洪保護區(qū)進洪量與蓄水量對比表

3.2.2流場分布 1)整體流場分布.通過對洪水淹沒的二維平面區(qū)域進行模擬，最終將模擬計算出的整體流場分布情況與DEM整體高程進行對比分析(圖8和圖9).對比圖8和9可知,淹沒區(qū)域流場分布均勻一致，而低洼地帶是較大流速主要集中的區(qū)域，根據(jù)洪水的自然流動原則，即從高到低流動，洪水呈現(xiàn)出的流動態(tài)勢較為準確.

圖8 和平村DEM地形示意圖

圖9 渝河遭遇100 a一遇洪水在和平村潰堤洪水最大流速示意

2)局部流場分布.通過對局部流場進行模擬分析，將模擬計算出的流場分布情況與線狀地物進行比較，對比分析結(jié)果見圖10和圖11(以和平村潰堤為例).結(jié)果表明，局部區(qū)域的流場分布均勻一致，且線狀地物的阻水效果比較明顯，洪水呈現(xiàn)出的流動態(tài)勢較為準確.

圖10 二維平面區(qū)域流場均勻分布

圖11 線狀地物阻水效果流場分布

3.2.3同一方案風險信息比較 如圖12～圖14所示，以和平村潰堤為例，將同一方案的DEM、洪水淹沒水深及洪水流速進行比較分析.結(jié)果表明，DEM地勢較高地區(qū)淹沒水深較小，相反，地勢較低洼地區(qū)淹沒水深較大，而對于處在地勢由高到低的區(qū)域的洪水，即DEM過渡地區(qū)的洪水，洪水流速偏大.這初步說明該方案計算具有合理性.

圖12 和平村潰堤DEM分布

圖13 和平村潰堤淹沒水深分布

圖14 和平村潰堤洪水流速分布

3.2.4不同方案風險信息比較 選定某一固定位置在同一時刻對不同方案的洪水流速、淹沒水深等洪水風險信息進行比較分析.以和平村潰堤為例，對渝河50 a和100 a一遇洪水進行風險信息比較.結(jié)果表明，同一位置同一時刻100 a一遇的洪水流速略大于50 a一遇洪水流速(圖15)；對于淹沒水深，100 a一遇洪水略大于50 a一遇洪水，比較結(jié)果見圖16.通過對不同方案風險信息進行驗證分析，可以發(fā)現(xiàn)本文模型計算較為合理.

圖15 同一位置同一時刻洪水流速比較

圖16 同一位置同一時刻淹沒水深比較

4 結(jié)語

綜上，對渝河(隆德段)洪水的耦合模擬，無論是模擬計算過程還是對結(jié)果的合理性分析，都對渝河區(qū)域洪水治理提供了一定的理論支撐，為水利防汛工作的開展提供依據(jù)，有利于防汛應(yīng)急預案的快速制定及部署.對于水利工程建設(shè)與管理相關(guān)部門來說，本文的模擬結(jié)果可以為渝河防洪治澇工程體系的規(guī)劃提供基本依據(jù)，便于開展有關(guān)洪水風險的宣傳及相關(guān)教育培訓等.