基于二維網(wǎng)格邊元設(shè)置河道方法的蓄滯洪區(qū)洪水演進(jìn)分析

喻海軍，吳俊峰，吳濱濱，孫宏圖，馬建明，高 佳

（1.中國(guó)水利水電科學(xué)研究院，北京 100038；2.水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038；3.河北省水利水電第二勘測(cè)設(shè)計(jì)研究院，河北 石家莊 050021；4.河北省防汛抗旱指揮部辦公室，河北 石家莊 050011）

1 研究背景

蓄滯洪區(qū)是我國(guó)江河防洪體系中的重要組成部分，是保障重點(diǎn)城市和地區(qū)防洪安全，減輕洪水災(zāi)害的有效措施[1]。蓄滯洪區(qū)洪水?dāng)?shù)值模擬對(duì)蓄滯洪區(qū)洪水管理、洪水損失評(píng)估以及流域洪水調(diào)度決策等都具有重要的意義[2-4]。近年來(lái)，國(guó)內(nèi)外學(xué)者在蓄滯洪區(qū)洪水?dāng)?shù)值模擬方面做了大量研究，取得了豐碩成果。王秀杰等[5]建立了河道與防洪保護(hù)區(qū)的全二維動(dòng)態(tài)耦合模型，實(shí)現(xiàn)了寬淺河道與防洪區(qū)的無(wú)縫銜接，較好地模擬了水流流態(tài)的變化。李大鳴等[6]采用一二維耦合水動(dòng)力學(xué)模型實(shí)現(xiàn)了大清河五洼滯洪區(qū)的實(shí)時(shí)洪水調(diào)度分析。Vorogushyn等[7]采用一維河道、潰堤和二維模型實(shí)現(xiàn)了潰堤洪水在蓄滯洪區(qū)內(nèi)的演進(jìn)模擬，對(duì)德國(guó)Elbe河蓄滯洪區(qū)內(nèi)的洪災(zāi)損失進(jìn)行了評(píng)估。一般來(lái)說(shuō)，大型蓄滯洪區(qū)的洪水模擬一般需要構(gòu)建一二維耦合模型[6-7]，但往往由于面積大、河流眾多等原因?qū)е履Ｐ蜁r(shí)效性比較差，給實(shí)時(shí)預(yù)報(bào)分析帶來(lái)一定的困難，因而十分有必要研究相對(duì)精準(zhǔn)且高效的數(shù)值模擬方法。本文采用基于二維網(wǎng)格邊元上設(shè)置河道的方法構(gòu)建一、二維耦合模型，并選擇海河流域最大的蓄滯洪區(qū)大陸澤寧晉泊蓄滯洪區(qū)為典型應(yīng)用區(qū)域，對(duì)本文建立的一、二維耦合模型進(jìn)行了驗(yàn)證與應(yīng)用，并分析了下滲因素對(duì)蓄滯洪區(qū)行洪的影響。

2 模型原理

2.1 一維模型河道洪水采用一維水動(dòng)力學(xué)模型模擬，控制方程為一維圣維南方程組：

式中：q為旁側(cè)入流；Q為流量；B為水面寬度；Z為水位；A為過(guò)水面積；α為動(dòng)量校正系數(shù)；Sf為摩阻比降。

采用基于有限體積法的Godunov格式對(duì)上述一維圣維南方程組進(jìn)行離散[8]。

2.2 二維模型地表洪水演進(jìn)采用二維水動(dòng)力學(xué)模型，控制方程采用守恒型的二維淺水方程：

式中：h、u、v、b分別為水深、x和y方向流速、底高程；Sb是底坡項(xiàng)；Sox、Soy分別是x和y方向的底坡比降；Sf是摩阻項(xiàng)；Sfx、Sfy分別是x和y方面的摩阻比降。

采用基于有限體積法的Godunov格式對(duì)上述二維淺水方程組進(jìn)行離散求解，其中Riemann問(wèn)題采用Roe格式的近似解進(jìn)行計(jì)算，底坡源項(xiàng)采用特征分級(jí)離散，保證模型的守恒性，阻力源項(xiàng)采用隱式離散提高模型的穩(wěn)定性，采用MUSCL空間重構(gòu)和兩步Runge-Kutta法使得模型具有時(shí)間和空間二階精度，所有變量都定義在單元中心[9-10]。

2.3 模型耦合一維河道和二維地表水動(dòng)力模型的水流交互采用邊元耦合的方式（如圖1所示），即在二維網(wǎng)格邊元上設(shè)置和建立一維河網(wǎng)模型，網(wǎng)格劃分時(shí)不考慮河道寬度，因而河道寬度不影響網(wǎng)格單元?jiǎng)澐?，可以極大的提高網(wǎng)格質(zhì)量，從而提高模型演算效率。相對(duì)于傳統(tǒng)分別構(gòu)建一二維模型，河道范圍內(nèi)不劃分網(wǎng)格或者全部采用二維模型計(jì)算，河道處進(jìn)行網(wǎng)格加密等方法，采用邊元設(shè)置河道來(lái)耦合的方式無(wú)論從網(wǎng)格劃分便捷性還是計(jì)算效率來(lái)看都得到了極大的提升，特別適用于在計(jì)算范圍內(nèi)存在較多對(duì)洪水演進(jìn)有影響的中小河道的情形。

圖1 二維網(wǎng)格邊元設(shè)置河道示意

圖2 河道與地面水流交互示意

如圖2所示，假設(shè)某時(shí)刻河道與地表網(wǎng)格水量通過(guò)側(cè)向交換的方式交換的流量為Ql，采用堰流公式近似計(jì)算交換流量的方法如下[11]：

其中：hmax和hmin分別采用下式計(jì)算：

式中：Zr、Zc分別為堰上、下游水位，分別取河道和二維網(wǎng)格單元的水位值；Ze為堰的高程，一般取河堤岸的高程；be為堰的寬度，一般取單元格與河道相連邊的邊長(zhǎng)。

3 模型構(gòu)建與應(yīng)用

3.1 區(qū)域概況本文選擇典型應(yīng)用區(qū)域大陸澤及寧晉泊蓄滯洪區(qū)位于河北省南部（如圖3所示），海河流域子牙河水系支流滏陽(yáng)河系中游，河北省邢臺(tái)市東北部，是全國(guó)第三大蓄滯洪區(qū)，也是海河流域第一大蓄滯洪區(qū)和關(guān)鍵防洪工程，總面積約2041 km2，總滯洪量35.52億m3，涉及邢臺(tái)市9個(gè)縣區(qū)，135.27萬(wàn)人。上游河流多處于太行山迎風(fēng)區(qū)，源短、坡陡、流急，洪水突發(fā)性、致災(zāi)性強(qiáng)。蓄滯洪區(qū)內(nèi)分布有眾多中小河流及堤壩，且堰閘、分區(qū)滯洪等調(diào)度規(guī)則復(fù)雜，導(dǎo)致一維河道與二維滯洪區(qū)洪水演進(jìn)過(guò)程復(fù)雜。大陸澤、寧晉泊蓄滯洪區(qū)按50年一遇洪水設(shè)計(jì)，從建國(guó)后至今60多年的蓄滯洪區(qū)運(yùn)用情況來(lái)看，大陸澤實(shí)際運(yùn)用4次，最高蓄洪水位33.41 m，最大蓄洪水量10.40億m3，蓄洪歷時(shí)57 d；寧晉泊實(shí)際運(yùn)用3次，最高蓄洪水位30.70 m，最大蓄洪水量29.27億m3，蓄洪歷時(shí)60 d。

3.2 模型構(gòu)建

（1）網(wǎng)格劃分。采用非結(jié)構(gòu)四邊形網(wǎng)格對(duì)整個(gè)研究區(qū)域進(jìn)行網(wǎng)格剖分，將區(qū)域內(nèi)高速公路、鐵路和河流作為內(nèi)部約束邊界，最終共劃分成約兩萬(wàn)個(gè)網(wǎng)格。在網(wǎng)格插值時(shí)，將堤防、高速公路、鐵路設(shè)置成線形阻水建筑物在網(wǎng)格邊元上予以考慮，將測(cè)量的實(shí)際高程作為阻水建筑物的高度，采用堰流公式計(jì)算建筑物兩側(cè)的過(guò)流水量。河道糙率取0.025，其它區(qū)域糙率取值范圍在0.040～0.050之間。基于網(wǎng)格邊元設(shè)置的河道分布如圖4（a）所示，網(wǎng)格局部放大如圖4（b）所示。

圖3 大陸澤寧晉泊蓄滯洪區(qū)位置示意

圖4 研究區(qū)域網(wǎng)格示意

（2）邊界條件。研究區(qū)域內(nèi)一共包含滏陽(yáng)河、留壘河、洺河、南澧河、順?biāo)樱ㄆ呃锖樱?、牛尾河、馬河（白馬河、小馬河和李陽(yáng)河）、泜河、午河（泲河）、北沙河（槐河）、小漳河和洨河等12條中小河道。在實(shí)際模擬計(jì)算時(shí)，河道上游入流采用流量邊界條件（典型洪水過(guò)程如圖5所示），下游艾辛莊樞紐的出流采用水位流量關(guān)系邊界（如圖6所示）。

（3）分區(qū)滯洪。大陸澤及寧晉泊蓄滯洪區(qū)被滏陽(yáng)河右堤和小漳河右堤劃成3個(gè)分區(qū)（如圖7所示），在實(shí)際調(diào)度運(yùn)行中，當(dāng)分區(qū)1的水位達(dá)到29.4 m時(shí)，由小南堤分洪口門(mén)向分區(qū)二分洪，當(dāng)分區(qū)二的水位超過(guò)29.4 m時(shí)由小漳河右堤分洪口門(mén)向分區(qū)三分洪。模型構(gòu)建時(shí)考慮了蓄滯洪區(qū)分區(qū)滯洪的調(diào)度規(guī)則，兩個(gè)分洪口門(mén)則采用潰口的方式進(jìn)行概化考慮，達(dá)到設(shè)定的條件（h＞29.4 m），自動(dòng)進(jìn)行分洪。

圖5 典型洪水過(guò)程（洺河20年和50年一遇設(shè)計(jì)洪水過(guò)程）

圖6 艾辛莊樞紐水位流量關(guān)系

圖7 分區(qū)滯洪示意

（4）模型驗(yàn)證。由于缺乏蓄滯洪區(qū)運(yùn)用的實(shí)測(cè)數(shù)據(jù)，無(wú)法對(duì)構(gòu)建的模型進(jìn)行較為準(zhǔn)確的率定驗(yàn)證，因而將模型計(jì)算結(jié)果與《子牙河系防洪規(guī)劃》（2008年）以及《大陸澤、寧晉泊蓄滯洪區(qū)建設(shè)與管理可行性研究報(bào)告》（2014年）中的模擬結(jié)果進(jìn)行比對(duì)分析，以對(duì)模型結(jié)果合理性進(jìn)行檢驗(yàn)?，F(xiàn)狀條件下，計(jì)算50年一遇洪水，采用位于蓄滯洪區(qū)上游的環(huán)水村（位于邢臺(tái)市任縣）以及下游艾辛莊兩處（位置示意見(jiàn)圖7）的最高洪水位進(jìn)行對(duì)比，結(jié)果如表1所示。總體上，本次構(gòu)建的模型與其它研究結(jié)果有一定差別，但水位差總體控制在0.15 m以內(nèi)，表明模型結(jié)果是可信的，差別主要來(lái)源可能是模型算法和概化方法不完全一樣。

3.3 模型應(yīng)用大陸澤及寧晉泊蓄滯洪區(qū)位于海河流域，屬于干旱半干旱地區(qū)，近幾十年由于地下水長(zhǎng)期處于超采狀態(tài)，地下水位持續(xù)下降，導(dǎo)致河道和地表行洪時(shí)下滲比較強(qiáng)烈，對(duì)洪水演進(jìn)（如洪峰、洪水淹沒(méi)時(shí)間等）有著重要的影響。由于以往針對(duì)大陸澤及寧晉泊蓄滯洪區(qū)的研究成果多沒(méi)有考慮下滲的影響[3]，從洪水調(diào)度管理的角度來(lái)看，不考慮下滲屬于偏安全的做法，存在一定的合理性，但由于蓄滯洪區(qū)一旦啟用，蓄水時(shí)間長(zhǎng)，下滲量就相對(duì)比較大，其影響不應(yīng)被忽略。本文著重對(duì)比分析了考慮和不考慮下滲情況下，蓄滯洪區(qū)洪水演進(jìn)情況的區(qū)別，下滲率取值參考程亮等[12]在海河流域的研究成果，穩(wěn)定下滲率取值12.3 mm/h。

圖8給出了考慮和不考慮下滲兩種情況下遭遇50年一遇洪水淹沒(méi)范圍。從圖8中可以看出，在t=50 h時(shí)上下游已經(jīng)有部分洪水從河道中溢出，但此時(shí)淹沒(méi)深度均較小，隨著時(shí)間增長(zhǎng)，淹沒(méi)范圍逐漸增大。總體上看，不考慮下滲時(shí)，50年一遇洪水需要啟用全部的三個(gè)分洪區(qū)，而考慮下滲時(shí)則只需啟用第一個(gè)分洪區(qū)，差別比較明顯。分析同一時(shí)刻考慮和不考慮下滲兩種情況，可以明顯看出同一時(shí)刻不考慮下滲淹沒(méi)范圍和淹沒(méi)深度均大于考慮下滲時(shí)，而考慮長(zhǎng)歷時(shí)洪水下滲對(duì)洪水淹沒(méi)范圍和淹沒(méi)深度都有較大影響。

表1 模擬結(jié)果對(duì)比

圖8 50年一遇不同時(shí)刻淹沒(méi)圖

圖9給出了20年一遇和50年一遇環(huán)水村附近洪水位變化過(guò)程。從圖9可以看出，考慮下滲時(shí)的洪峰值明顯小于不考慮下滲時(shí)，數(shù)值上分別減少了0.42 m和0.41 m，另外，考慮下滲時(shí)的退水時(shí)間與不考慮下滲時(shí)相比分別減少了195 h和296 h。從整個(gè)蓄滯洪區(qū)來(lái)看，以遭遇50年一遇洪水為例，不考慮下滲時(shí)蓄滯洪區(qū)需要約60 d的時(shí)間將洪水基本排盡，考慮下滲時(shí)，則只需要約36 d時(shí)間，減少了近24 d時(shí)間。

圖9 環(huán)水村洪水變化過(guò)程

4 結(jié)論

（1）基于在二維邊元設(shè)置河道方法，建立了大陸澤及寧晉泊蓄滯洪區(qū)一二維耦合的地表和河道洪水?dāng)?shù)值模型，在與其它模型計(jì)算結(jié)果進(jìn)行對(duì)比驗(yàn)證時(shí)，水位差總體控制在0.15 m以內(nèi)，表明該模型結(jié)果是可信的，在此基礎(chǔ)上對(duì)20年一遇和50年一遇洪水進(jìn)行了模擬與分析。應(yīng)用結(jié)果表明基于二維邊元設(shè)置河道的一二維耦合模型能夠很好的應(yīng)用于中小河流較多的蓄滯洪區(qū)，具有廣闊的應(yīng)用前景。

（2）對(duì)比分析地表下滲對(duì)蓄滯洪區(qū)洪水的影響可見(jiàn)：考慮下滲作用時(shí)，蓄滯洪區(qū)50年一遇洪水只需要啟用1個(gè)分區(qū)，而不考慮時(shí)則需要啟用3個(gè)；考慮下滲時(shí)退水時(shí)間與不考慮下滲時(shí)相比也相應(yīng)減少了近24 d。表明干旱半干旱地區(qū)下滲對(duì)蓄滯洪區(qū)行洪具有較大的影響，在蓄滯洪區(qū)實(shí)際調(diào)度應(yīng)用時(shí)，應(yīng)該考慮下滲因素對(duì)洪量減少的影響，以較好發(fā)揮蓄滯洪區(qū)的效用。